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ABSTRACT 


A computer program has been developed to solve the low speed flow 
around airfoils with highly separated flow. A new flow model, which 
was suggested by Zumwalt , included all of the major physical features 
in the separated region. It was suggested by experimental airfoil 
studies made in the WSU low speed wind tunnel. Flow visualization tests 
also were made which gave substantiation to the validity of the model. 

The computation involves the matching of the potential flow, bound- 
ary layer and flows in the separated regions. 

The potential flow program was available from the McDonnel 1 -Douglas 
company. Head's entrainment theory was used for boundary layer calcula- 
tions and Korst's jet mixing analysis was used in the separated regions. 
A free stagnation point aft of the airfoil and a standing vortex in the 
separated region were modelled and computed. The separation location 
and pressure were found iteratively without a priori specification. 

A GA(W)-1, 17% thick airfoil, at three angles of attack and two 
Reynold's numbers, was used for the analysis since experimental data 
were available. The surface pressures resulting from the computation 
matched very well with experimental data. In particular, separation lo- 
cations and pressures were nearly identical with experimental values. 
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I. INTRODUCTION 


Separated flow behind wings has been of interest to researchers from 
the time it was first discovered that it causes stalling of aircraft. 

The complex nature of the problem hindered the analysts in obtaining a 
solution to the problem of separation because solving the Navier Stokes 
Equation was an insurmountable task before high speed computers were 
available. As technology improved, making possible short take offs and 
landings, separated flow research gained more importance. Research in 
the past has been mainly experimental due to the complexity involved in 
theoretical analysis. However, with the advent of modern computers it 
has been possible in principle to solve Navier Stokes Equations, but 
with the available memory and speed the cost of analysis has been pro- 
hibitive. Several mathematical models have been proposed and solved 
by computers which did not involve a complete solution of the Navier 
Stokes Equation, but used simplified boundary layer methods for the 
viscous flow, some empirical relations, and a few assumptions. Most of 
the models hitherto proposed have considered the separated region either 
as extending to infinity or as a bulbous region. The specific details 
in the separated region have not been considered.' Axisymmetric separ- 
ated flow analysis is available but is not applicable to the highly un- 
symmetric separated region behind an airfoil at high angles of attack. 

The present separated flow model was suggested by Zumwalt (Ref. 1) 
based on experimental measurements of separated flow on a GA(W)-1 air- 
foil in the WSU 7' x 10' subsonic wind tunnel. It takes into account 
the details in the separated region. The velocity field pattern pre- 
dicted by this model closely matched the experimental measurements thus 


1 



providing a basis for the validity of the model. Flow visualization 
studies of the separated flow behind an airfoil gave qualitative sub- 
stantiation. 

The analysis involves a computer solution for low velocity flow 
around and behind an airfoil with massive separation. It required 
first an invisicid flow analysis, and second a matching with all pos- 
sible viscous interactions. 

Typical of successful computation programs for attached flow on 
airfoils are References 6 and 7. Ref. 6 is the direct method where 
specified geometry produces pressure distribution, and Ref. 7 is the 
inverse method where specified pressure distribution gives geometry. The 
McDonnel 1 -Douglas Mixed Boundary Condition (Ref. 3) program will permit 
one to either (a) supply the surface profile and obtain the adjacent 
flow conditions or (b) supply the pressure of a streamline and obtain 
the streamline location. This provides a single method for treating 
both attached and separated flows. 

Separated flow occurs when the flow leaves the surface of the air- 
foil due to an adverse pressure gradient. The location of the separa- 
tion point plays the important role of dividing the flow regimes. It 
is a point of zero surface shear and is dependent on the pressure gra- 
dient, angle of attack and the nature of the boundary layer. 

Several investigators have developed criteria to predict the sep- 
aration point using analytical and empirical methods. Reference (4) 
has described some of these methods and evaluated the methods by com- 
parison with experiment. There is no agreement among the different 
methods and all of them draw from experimental data in order to be able 
to predict the separation point. Predictions were also dependent on 


2 



the method used for boundary layer calculations and the method of usinq 
experimental data for the empirical evaluation. Thus it is clear that 
there is still some confusion in the empirical methods for predictinq 
separation. 

The shape factor H (H = 6*/e) has been used as a guide for deter- 
mining the separation point. Earlier investigators have successfully 
used this method for prescribing separation, and in the present analy- 
sis, separation is specified by prescribing a value of the shape fac- 
tor, H. 

The existence of two kinds of flow, namely laminar and turbulent, 
necessitates the distinction between laminar separation and turbulent 
separation. Airfoil separation has been classified into three cate- 
gories by McCullough and Gault (Ref. 5). 

(1) Tra iling edge separation . This is essentially a turbulent 
separation at the trailing edge moving upstream with angle of attcick 
increase. The turbulent separation can either be preceded by transi- 
tion from laminar to turbulent or by a laminar separation and with 
turbulent reattachment. The latter is the case of a short bubble with 
transition. 

(2) Leading ed ge sepa ration: laminar flow separation near the 

leading edge without any reattachment. 

(3) Thin airfoil separation : laminar flow separation near tne 

leading edge with flow reattachment (long bubble) at a point which move 
downstream with increase of angle of attack. 

The present analysis assumes that the flow is turbulent very neat’ 
the leading edge and thus considers only t urbulen t sepa ration . Laminar 
separation is not included since it makes the analysis more complicated 
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and is unlikely in practical low speed aircraft wings. Also, the main 
object was to analyze the flow model with attention to the details in 
the separated region. The flow is assumed to be steady, incompressible 
and two-dimensional. 

Head's entrainment method was adopted to calculate the turbulent 
boundary layer characteristics and the separation point and was used in 
an iterative interaction with the inviscid flow solution, until conver- 
gence was achieved. 

Korst's separated flow analysis was used to map the flow from sep- 
aration point to the trailing edge station of the airfoil for the upper 
surface. Mass flows into and out of the separated region behind the 
airfoil is also determined by the above analysis and from boundary 
layer theory. 

The two separated flows, namely the top surface separated flow and 
the bottom surface separated flow, and the free stagnation point at 
which they meet formed the separation bubble behind the airfoil. 

The free stagnation point behind the airfoil was located at a dis- 
tance of a third of the distance between the separation point and the 
trailing edge from the trailing edge. This was an assumption based on 
the WSU experimental measurements (Ref. 2). The vertical position of 
the free stagnation point was not fixed. 

The details of flow in the separation bubble circulatory flow was 
determined by velocity distributions of Korst's model and an assumed 
profile for the reversed flow based on experimental data. 

A computer program has been developed to solve the separated flow 
model on a two dimensional airfoil for incompressible subsonic flow. It 
determines the pressure distributions on the airfoil and maps the separated 
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flow region. It uses the McDonnell -Douglas potential flow program for 
the inviscid analysis. Head's entrainment method and Korst's separated 
flow model for the viscous flow and the separated region respectively. 

The conditions to be satisfied for a valid solution are discussed in 
detail in Chapter III. 

Flow visualization experiments were made for a qualitative obser- 
vation of the flow details in the separated region and they generally 
indicated that the assumptions in the flow model were realistic. 

The computer program was used to calculate the pressure distribu- 
tion around a 17% thick GA(W)-1 wing (Fig. 1) at three angles of attack. 
The results were compared to experimental data to verify the model and 
computational method. 

This research was conducted under a grant from NASA Langley Research 
Center; Grant No. NSG 1192. 
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II. PREVIOUS RELATED WORK 


The analysis of separated flow on airfoils has been approached in 
two ways: one, by a numerical solution of the complete Navier Stokes 

equations, and the other by assuming a physical model for the separated 
region and then solving for the mathematical solution. 

Numerical solution of Navier Stokes equations for separated flow 
has been attempted by Thames et al (Refs. 8) on arbitrary airfoils. The 
analysis reported is only for laminar flow and for very low Reynolds num- 
ber. Fig. 2 shows a typical streamline pattern obtained by this method 
of separated flow on an airfoil. Higher Reynolds number analysis will 
require prohibitive length of computer times. Turbulent modelling in 
separated regions is the most important aspect of separated flow and 
this is not achieved in this method. 

Jacob's (Refs. 9,10,12) idea (Fig. 3) of using a source or source 
distribution in the aft region of the airfoil to form the separated re- 
gion has been adopted by a number of researchers with modifications. 

Hahn et al, Bhately and Mcwhirter, Farn et al, are some who have reported 
analyses based on this idea (Refs. 4, 11,13). The main differences in 
all these models are in their methods of finding the source distributions 
which satisfies the boundary conditions and of choosing the location of 
the separation point and pressure. 

In all these models the separated region is considered to extend to 
infinity which is quite contrary to experimental observations. Jacob 
(1975) proposed a model recently which closes the region by using a 
sink at a point downstream of the trailing edge (see Fig. 6). The de- 
tails and some of the results are shown in Figs. 4 and 5. 
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Extension of a conventional boundary layer analysis into the 
separated region has also been attempted (Refs. 14,15) for the repre- 
sentation of separated flow. But, as one can see, the credibility of 
the boundary layer assumptions are lost when the boundary layer thickness 
becomes very large. 

In addition, the available methods for separated flow analyses 
have not considered the effect of reverse flow on the pressure distri- 
bution of the airfoil. The separation pressure or its location has 
been assumed. 

Based on experimental investigations (Figs. 7,8,9,10) in the 7' x 10' 
WSU low speed wind tunnel (Ref. 2), a model will now be presented which 
takes into account more details in the separated region. One of the key 
factors not considered in the earlier models is the mass recirculated in 
the separation bubble. 
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III. DESCRIPTION OF THE FLOW MODEL 


A sketch of the flow model is shown in Fig. 11. The. flow has sep- 
arated from the upper surface and leaves the lower surface at the trail- 
ing edge. The jet mixing sheets starting from the two airfoil separa- 
tion points coalesce to form the separation bubble behind the airfoil. 
These two jets entrain air from the dead air region (near-wake). The 
entrained air has to be replaced by backflow of air which must be sup- 
plied from somewhere. If S and S' (Fig. 11) are the two separation 
streamlines, it can be seen that the amount and width of flow entrained 
are growing in the downstream direction, consequently, the space for 
backflow decreases and the demand for it increases. Since this cannot 
continue very far, a termination is required of the near wake recircu- 
lation region, and a stagnation point is formed. At this point the 
two streams rejoin defining the end of the separation bubble. This 
near wake free stagnation point is termed a "saddle point" since pres- 
sure distribution resembles a saddle. 

The entrained masses of the two jets are not the same since their 
lengths and velocities are different. Generally, the upper one will 
entrain a larger mass. Therefore the stagnating streamlines are not 
necessarily S and S' at the saddle point. 

Details of the separated regions are shown in Fig. 12. it shows 
two other streamlines R and R' stagnating at the saddle point provid- 
ing passages for the flow to enter and leave the separation bubble. 

The mass flow through the corridor between R and S should be the same 
as between R' and S'. Further, R and R' do not terminate at the saddle 
point but must continue upstream and downstream, the required mass 
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conservation in the region requires the formation of two bubbles and 
an S-shaped corridor. Based on experimental data, constant pressure 
is assumed in the whole region except in the neighborhood of the sad- 
dle point, where a high pressure wedge will form along the RR' line 
and extend into the recirculating regions, turning back the low veloc- 
ity flows. 

Analysis of flow in the separated region would require a nearly 
constant pressure along S and S' in order to be able to apply turbu- 
lent jet mixing analyses for constant pressure regions. Experimental 
data have indicated that this is true on the upper surface but not in 
the vicinity of the saddle point. 

Viscous effects can be ignored in the neighborhood of the saddle 
point and all velocity changes considered as being due to the pressure 
gradients. Thus the region is divided into viscous dominated and pres- 
sure dominated regions. This follows the classical approach originated 
by Prandtl for boundary layers and more recently applied successfully 
to separated and reattaching flows in the Chapman-Korst (Ref. 12) mix- 
ing models. 

The jet mixing theory is required for computation of turbulent mix- 
ing entrainment rates. The Gbertler exchange coefficient model, as 
adapted by Korst, was considered to be best for a first attempt due to 
its successful application to other plane flow problems and the avail- 
ability of mass and momentum integral tables for these flows. 

The trailing edge plane is assumed to divide the constant pressure 
region of the separated flow and the region of pressure rise to the 
near stagnation point. This plane is also the location for satisfying 
the mass continuity. 
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Fig. 18 shows a schematic diagram of the velocity profile at the 
trailing edge of the airfoil. The velocity profile consists of three 
segments: 

(1) The error function profile of the upper high velocity flow. 

(2) A constant velocity reverse flow region. 

(3) A third degree parabola to join (1) and (2). 

The Korst profile is adopted until the point where the velocity is half 
the value at the outer edge of the shear flow. The parabola matches 
slopes and velocity at the other two flow segments. The matching lo- 
cation for the reversed flow is chosen from experimental data, as will 
be explained in Chapter V. 

The complete solution of a wing requires the mating of the separated 
region to the potential flow and boundary layer flow. 
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IV. EXPLORATORY EXPERIMENTS 


A. ELECTRICAL ANALOGY 


An electric conducting table analogy experiment was attempted ini- 
tially in the hope that it would indicate the conditions required for 
the formation of the separated region with a potential flow model and 
that it would also aid parameter selection for the computer program. 

A GA(W)-1, 13% thick, 10 inch chord airfoil was used for this purpose. 

Fig. 13 shows the apparatus. It consists of the conducting paper laid 
flat on a table with the silver-paint airfoil in the center of the 
paper. The ends of the paper are firmly held by conducting rods or 
angle sections. The airfoil is oriented on the paper such that the 
streamlines are parallel to the conducting rods. Electrical leads 
buried in the paint are connected to a potentiometer to vary the volt- 
age of the airfoil . 

The model simulation was attempted by placing circular brass discs 
of h, inch diameter at positions where the vortices were predicted. The 
voltages in the airfoil and the two discs were varied to represent vari- 
ous values of circulation. The flow was simulated by reducing the volt- 
age of the airfoil so that the trailing edge stagnation streamline was 
displaced and formed on the upper surface of the airfoil. This corres- 
ponds to the reduction of circulation on the airfoil. The voltage of 
the disc near the airfoil was then adjusted so that the streamlines very 
near the lower surface were diverted to form the S-shape as represented 
in the model. Fig. 14 shows the streamline shapes resulting in the 
electric analog. This indicated that vortices can be used in the computer 



program to form the separated region but the separation cannot be simu- 
lated satisfactorily for the potential flow model, but rather a boundary 
layer must be included in the attached flow program. The usable result 
from this experiment turned out to be the fact that small movements of 
the downstream saddle point did not change the flow pattern appreciably. 
This was an important assumption made in the computer program since the 
vertical position of the saddle point is not specified. The idea of 
using vortices to form the separated region was abandoned after it was 
known that the potential flow program to be used in the analysis could 
provide the displacement surface of the separated flow. 

B. FLOW VISUALIZATION STUDIES 

Flow visualization experiments were conducted in a small, 6" x 14", 
low speed wind tunnel at a Reynolds number of 0.3 x 10^. A 10 inch 
chord, two-dimensional, GA(W)-1 wing was used. The wing was held at 
a fixed angle of attack, supported by the sides of the tunnel. A thin 
aluminum plate was mounted vertically on the wing at midspan so that it 
was parallel to the flow. The plate was smeared with a mixture of lamp- 
black and kerosene. A few flow photographs are shown in Figs. 15 and 
16. The main observations are as follows: 

(1) The wake closes behind the airfoil to form a bubble shaped 
region with the free stagnation point very close to the trailing edge 
as the earlier measurements of Ref. 2 had indicated. 

(2) The recirculating flow in the near wake forms a fairly large 
unsymmetric vortex which is clearly seen in the photographs. 

(3) There is an upward flow from the lower wake of the airfoil 
flowing upstream in separation bubble and turning back to join the 


12 



main stream. This S-shaped, lower- to-upper flow may not be clearly 
visible in the photographs, but was easily detected during the tests. 

The flow visualization tests confirmed that the assumed features 
in the model were present and thus substantiated the validity of the 
present model . 
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V. DESCRIPTION OF ANALYSIS OF THE FLOW 


The analysis of separated flow over an airfoil has been divided 
into two sections, an outside problem and an inside problem, each of 
which is solved separately, then matched for a complete solution. Fig. 
17 shows the geometry of the outside problem. Velocities and pressures 
along the attached-flow surfaces of the airfoil were calculated by a 
computer program adapted from the McDonnel 1 -Douglas Mixed Boundary Con- 
dition program. A boundary layer computation was made along the sur- 
face to indicate the separation point and provide the displacement thick 
ness. These were iterated until the solutions became stable. The pres- 
sure of the separated region was assumed to be equal to separation point 
pressure from the separation point to the airfoil trailing edge. A para 
bolic increase to the free stagnation point pressure from the trailing 
edge pressure was assumed. 

The separation point location was determined by the boundary layer 
routine. Fig. 18 shows the inside problem. The Korst jet mixing analy- 
sis was used to determine the mass entrained from the separated region 
by the upper surface jet sheet. The two stagnating streamlines, R' and 
R, and the rear free stagnation point pressure were determined from a 
balance of the mass inflow and mass outflow of the separation bubble. 

The mass inflow into the separation bubble from the lower surface of 
the airfoil was calculated by using a power law velocity profile for 
the turbulent boundary layer at the trailing edge station. 

A. POTENTIAL FLOW SOLUTION 

The Mixed Boundary Condition flow program of McDonnel 1 -Douglas 
Company (MCAIR) (Ref. 3) is a modification of the wing body analysis 
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program developed by Woodward (Ref. 18). The configuration of the air- 
foil is divided into a number of panels on the chordline. The effects 
of thickness, camber and the angle of attack are represented by planar 
source and vortex singularities. The boundary condition, the Neuman 
and the Kutta conditions determine the strengths of the source and vor- 
tex singularities. They form a system of linear equations which are 
then solved for the singularity strengths. The program allows the spe- 
cified boundary conditions to be given either as surface geometry or 
surface pressure distributions. The equations are solved by a routine 
using Gauss elimination to obtain the pressure distributions and sur- 
face configuration of the airfoil and the separation bubble.. The pre- 
sent analyses are made on a GA(W)-1 general aviation airfoil (Ref. 19). 

The geometry conditions required by the program are the slope of 
the camberline and the slope of the airfoil thickness distribution. In 
the present case, the airfoil geometry was known. In a more general 
case only the airfoil thickness distribution and the camberline are 
specified at a number of points on the chordline. A subroutine in the 
program prepares the geometric boundary conditions in the required form 
and at required stages of the program. The Neuman condition requires 
that the velocity be tangent to the surface, implying no flow across 
the physical boundaries. The Kutta condition determines the unique cir- 
culation around the airfoil. This is satisfied by specifying upper and 
lower pressures to be equal at the trailing edge. This enables Kutta 
condition to be satisfied for airfoils with blunt trailing edges as in 
the present case. However, when separation is present, the point at 
which the Kutta condition is to be applied is generally not clear. In 
the present case it is satisfied by specifying the same pressures for 
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the separation point and the lower surface trailing edge point. The 
points at which the singularities are applied are important since when 
discrete singularities are applied on finite sized panels there is a 
mathematical singularity at each edge of the panels and the velocity 
or pressure calculated is erroneous. This problem is avoided by apply- 
ing the boundary conditions for the vortex singularity at an interme- 
diate point and choosing this "control point" in such a way that the 
resulting solution at this point is as near as possible to the correct 
one. The optimum control point in the present case is at 85% of the 
panel length, as suggested in Ref. 3. The accuracy of the overall 
solution will depend upon the size and number of panels. An improve- 
ment is seen if the control point corresponds as closely to the trail- 
ing edge as possible. But on account of numerical instability there 
can be no sharp disparity in adjacent panels. Hence, a trade-off is 
established between the accuracy and the cost of calculation. Simi- 
larly, the nose region should also be represented by a larger number 
of smaller panels on account of its high curvature. The panel length, 
in the leading and trailing edge portions were chosen to be 1% of 
the chord, and the lengths increased to 5% in the center of the airfoil. 
The panels in the near wake were also 1% long. 

B. BOUNDARY LAYER ANALYSIS 

Viscous flow interactions are introduced by adding the displace- 
ment thickness of the boundary layer to the original airfoil shape. 

The new pressure distribution on the augmented airfoil is then deter- 
mined by the potential flow program. The boundary layer characteristics 
of the augmented airfoil give a new value of displacement thickness. 


16 



This is again added to the original airfoil and the iterative process 
is continued until the pressure distribution settles down to within 0.01. 

When separation occurs on the top surface, the iterations are con- 
tinued only up to the separation point. Since the potential flow pro- 
gram also provides the separation streamline, the displacement thickness 
at the separation point is added to the ordinates of the separation 
streamline up to the trailing edge to obtain the displacement surface. 

A parabolic pressure distribution was assumed for the region after the 
trailing edge up to the rear stagnation point for both the upper and 
lower surfaces. To start the iteration, a value for the rear stagna- 
tion pressure was assumed. 

Since the objective of the project was to show the feasibility of 
the mathematical model, sophisticated boundary layer analysis .was not 
used. Instead, the flow was assumed to be turbulent from the leading 
edge. The momentum integral method was used because of its simplicity 
and adaptability to iterative calculations. 

Head's entrainment method (Ref. 20) of calculating the turbulent 
boundary layer characteristics was chosen as being sufficiently accur- 
ate without undue complexity. 

The momentum integral equation for incompressible, two-dimensional 
flow in the integral form is: 


de _ ^f 
dx 2 


ue 


dug 

dx 


(H + 2) 


2.1 


where 



0 = Momentum thickness of the boundary layer 
Cf = Skin friction coefficient 
Ug = Velocity at the edge of the boundary layer 
H = Shape factor = 6*/9 
6* = Displacement thickness 


Head introduced the concept of the mass entrainment to the boundary 
layer. He argued that the rate of change of mass within the boundary 
layer was a unique function of the velocity defect. He derived a method 
for calculating simultaneously the development of the momentum thick- 
ness 0 and a quantity a which is referred to as the mass flow thick- 
ness. 



2.2 


An auxiliary equation is obtained by considering the rate at which the 
turbulent boundary layer entrains fluid from the free stream. 

dA _ p A ^^e 

dx u^ dx 2.3 

The non-dimensional entrainment parameter F is a unique function of 
another shape factor. Hi = a/0. Head obtained empirical relationships 
for F and relations between H] and the familiar shape factor H = S*/d. 

F = 0.0306 (Hi - 3.0)'°-®^^ 2.4 

2.5 
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Green (Ref. 20) obtained a linear relation for F by cross plotting in 
terms of H. 


F = 0.025H - 0.022 


2.6 


He also considered a relation between H and H] which was in better 
agreement with experiments. 


= 3.3 + 


0.9 


(H-l)‘*/3 


2.7 


The auxiliary equation is simplified by writing equation 2.3 as 


u do _ dHi 

^1 d7 " ® IT ■ ■ ^1 


0 dUg 

dx 


Using the relation between H and (Eq. 2.5), the auxiliary equation 
becomes: 


e^- H 


(h2- 

Ug dx 


H-1 

2 


|]h- 1) F - HCf] 


2.8 


The momentum integral equation, together with the auxiliary equation 
and relations for the local skin friction provide a step-by-step method 
for calculating the development of an incompressible turbulent bound- 
ary layer. Two skin friction relations were considered, those of 
Felsch and Ludwig-Tillman. The two expressions are: 
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Cp = 0.246 ^g-0.678H Ludwig -Tillman 

Cp = 0.058 (0.93 - 1.95 log Re"°-^^®... Felsch 

Both have wide acceptance and when used in the present program gave 
very similar results. However, the Ludwig-Til Iman expression seemed 
to amplify a numerical instability tendency at one point in the compu- 
tation development, while the Felsch did not. Hence, the Felsch ex- 
pression was retained. 

C. SEPARATED FLOW 

Based on the experience of other investigators, separation was 
assumed to have occurred when the value of H reached about 2 on the 
upper surface of the airfoil. Even though H values as large as 2.6 
have been measured experimentally, the higher values of H and the rapid 
boundary layer growth produced large induced slopes at the panels caus- 
ing severe instabilities in previous analyses (see Ref. 21). Although 
specification of an exact value of H for separation is not possible, 
most of the integral methods have assumed H from 1.8 to 2.4. Even 
though the range of H seems large, the separation location does not 
vary as much since close to separation the shape factor increases quickly. 

1 . Jet Mixing Analysis 

Separation is assumed at the end of any panel in which it occurs. 
Korst's theory, developed for a constant pressure, turbulent mixing of 
an isoenergetic free jet, was adapted for the incompressible case here 
to model the flow on the upper surface from the separation point to 
the free stagnation point. 
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Korst assumed similar velocity profiles in the free shear layers 
for fully developed flows and proposed that the velocity profile could 
be modelled by 


4, = 1 

Ue 2 


1 + erf (o|-) 


^(1 


erfn) 


where a is the jet spreading rate parameter. A value of a = 12 is well 
established for subsonic flow. For the viscous shear layer an intrinsic 
system of coordinates defined by the center of the mixing region (i.e., 
y = 0 at u = %ue) creates a shift y^, between intrinsic and inviscid 
coordinate systems. This shift is determined by the use of the stream- 
wise momentum equation. 

Equating the momenta below the e streamline (Fig. 21) for the two 
sections gives the shift ym. 


ye 



= Pe ^ (1 - cj)l2^ 


where ye = edge of the shear layer, and n = a 


X ’ 



00 
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A jet boundary streamline (Fig. 19 ), S, which separates the mass origin- 
ally flowing at the separation point from that entrained by jet mixing 
from the dead air region is found by equating the mass for the given 
velocity profile with that of undisturbed flow. 

ye 

M - Pe Ug (yg - y^j,) = Ug2 — (og - = j pudy 

ys 


-Pe Ug 0 (1 - Ce^ ) (Ilg - 


ne 


where 


I 


(|)dn 


ns 


A - 


2^2 


and 


/ 


ij>dn 

A - ce^<J>^ 


The mass flow integrals are tabulated for various values of ^ 

(<j) = u/ue) and ce (Crocco Number). Thus the mass flow between any two 
streamlines in the mixing region could be obtained by taking the dif- 
ference between two integrals corresponding to them and multiplying it 
by the appropriate variables of the flow. 

The boundary layer at separation can be replaced by a jet mixing 
profile having the similar characteristics. This is done by locating 
a virtual origin for the mixing which gives the resulting jet mixing 
profile which is the same as those of the actual boundary layer. The 
virtual origin is displaced upstream of the separation point by Xq 
and yo on the intrinsic coordinate system. This method was developed 
by Hill (Ref. 23). The details are shown in Fig. 19. 
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b: the actual streamline corresponding to the displacement 

surface from the separation point. 

X: inviscid streamline for the Korst flow starting from 

(xo>yo)* Note Xand b are parallel. 

x: intrinsic coordinate axis corresponding to a velocity of 

half the value in the inviscid stream adjacent to the dissipative 
region 

S: streamline which separates the mass originally flowing at 

the separation point (or more precisely at the virtual origin) from 
that entrained from the dead air region. 

The expressions for Xq and yg are 


xo 


(1 - c2e)I], 


From (Ref. 23) 


% 


yo = e + 6* 


where 6* and 0 are the displacement thickness and momentum thickness 
at the separation point* 


and 
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or incompressible flow: 


Cg^ -*■ 0. 


Xo = 



ns ns 

Ils = / ‘I'Cln = f f + erfn)dn 

->CO —CO 

4)5 = 0.61632 for Cg2 = 0 From (Ref. 24) 

From the tables of I-| integrals we obtain: 

II 5 = 0.399 
12 

Xq “ Q 39g ® ~ 

and 

yo = e + 6* 

2. Stagnation Streamline Determination 

Initially a value for the rear stagnation pressure is assumed, 
based on experimental data, to start the calculations. This enables 
the determination of the two displacement surfaces from the upper and 
lower separation points. There is an upper limit to the choice of the 
reattachment pressure since it is the stagnation point for the flow 
inside the separation bubble. Referring to Fig. 19, the S streamline 
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which separates the primary mass from that entrained from the separated 
region will have the highest possible stagnation pressure in the separa- 
tion bubble. 


■ = ~ . Ps P^s^ - P- 

S Q oo Q oo 

^ PsEP ' Pg. ^s^ ^e^ 

qo= uj- uj. 

2 

^Ps " ^PSEP uf 



CO 


Alternately, the cfioice can be based on the experimental value from Ref. 
14. The initial value of Cpp was taken as 0.0 to begin the calculations. 

The mass leaving the separation bubble M|j is obtained by consider- 
ing the flow between the S streamline and the R streamline, which stag- 
nates at the stagnation point. The mass entering, M|_, is obtained by 
considering the flow at the trailing edge between the lower surface of 
the airfoil and the stagnating streamline R', see Fig. 18. 

From Korst's analysis: 


U.Xti: 

= ^(11, - Hr) 


= Ue (100 - xsEP + Xq) (n. - Hn)/a 


and lit- can be found from the tables for known values of and 
Since the proper <j)p is determined by iterating the mass flows entering 
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and leaving the separation bubble, a parabolic fit for I-]j^ as a function 
of for a working range of n values was determined from the mass flow 
tables. The expression is: 

= 0.11402 - 0.20457 + 0.0817 <j>^^ 

Values calculated by this expression matched within 0.03% for the range 

n = 0 to n = n^* 

Ml is calculated from boundary theory, as follows. Here, the prime 
indicates lower side conditions. 



Assuming a power law for the boundary layer profile: 




6 1 
= / (ue - u')dz’ = ue6‘ / 
0 0 


(1 - 4>')d 





★ 

6 ' 


JL.) = _L 

n+1 / n+1 


Now 


6 

UeV = y* u'(ue - u')dz = ue^6' 
0 



(i>')2]d 



0' _ n 
6' (n+1 )(n+2) 
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A relationship between the shape factor H' and n can be determined: 




H’ = 


n+1 


n+2 


0* 


(n+D(n+2) 


n = 


H' 1 

For known values of H', n can be found. 


Hence Ml is calculated as: 


/ ,\l/n 




" K/6 / V 

/ (^) 


^ d 


5' = 6'*(n + 1) 


z ’ D I / p* ■ 


R76^ 


Ue(n + 1)<S f 
0 


1/n 


(n+1 )6^ 


,(n+l)6' 


p 




where 


n = 


H’ - 1 


Since the streamlines R and R’ stagnate at the same point from 
the same static pressure, their velocities must be equal: 

Up = K 

The value of u^ is iterated until My = Ml- 
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3. Recirculating Mass Balance 

The mass balance for the separation bubble is calculated at the 
99.85% chord station, which is the control point of the last airfoil 
panel. The velocity profile is made up of four segments. Figure 24: 

a) From the R streamline to the u=0.5 Ue streamline, the error 
function jet mixing profile is used. Thus the upper half of 
the ordinary free jet profile is retained. 

b) Experiments show a constant-velocity reverse flow, u^^ (i.e.> 
forward on the airfoil) region. This is assumed to exist from 
the inside of the error function profile ( defined as u=0.01 
or n = -2) to the augumented airfoil surface. The value of this 
reverse flow velocity will be discussed later. 

c) Between (a) and (b), a third degree parabola is placed which 
matches values of slopes of both the (a) and (b) profiles. 

d) The airfoil upper surface is augumented by the displacement 
thickness of the lower surface trailing edge. This pictures 
the boundary layer as swirling almost unchanged around a small 
separated bubble at the trailing edge. 

Since the evidence for a constant-velocity reverse flow profile is great, 
no logical and simple model can be suggested to give a"core" flow from 
reversal of shear flows. A purely empirical choice was accepted as nece- 
ssary and a value of Ur=0.2 Ue, was derived from examination of several 
GA(W) wing flow measurements. The sensitivity of the results to this 
choice will be discussed later. 

M - ^SEP^TE , 

'^STREAMWISE ' a " ^n=0^ 

n 

^^REVERSE = r udz + • Q • u^ep 

j ^SEP 

n — 2 
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where no = value of n for u = 0. 


The mass balance in the recirculating region was achieved by chang- 
ing the value of the shape .factor H for separation and iterating until 
convergence is reached. The pressure distribution corresponding to this 
condition was accepted as the solution of the analysis. 

D. ASSUMPTIONS IN THE PRESENT ANALYSIS 

1. The flow is assumed to be incompressible and steady 

2. The boundary layer is fully turbulent. 

3. The rear stagnation point location is assumed at a chordwise dist- 
ance of one-third the distance from the separation point on the 
upper surface to the trailing edge. This was based on the experi- 
mental results of Ref. 2. 

4. The rear stagnation point pressure coefficient value assumed to 
start the calculations. This was also based on the experimental 
data of Ref. 2. This value is, of course, replaced to form 
convergence 2. 

5. A constant 'core' velocity of the reverse flow in the separation 
bubble equal to a value of 0.2 times that of the velocity at the 
edge of the shear layer. 


29 



VI. DESCRIPTION OF THE COMPUTER PROGRAM 

The program is titled "SEPFLO" and is coded in fortran language 
to operate on the IBM 360 or other compatable models. It is divided 
into three main sections, namely, the potential flow, the boundary 
layer and the separated region. 

The potential flow part was adapted from the McDonnell Douglas 
Aircraft Company's Mixed Boundary Condition Program (Ref. 3) developed 
from the earlier program of Reference 18. The McDonnell Douglas prog- 
gram is still a proprietary item and hence it will not be discussed in 
detail here. Henceforth it will be referred to as the "Potential flow 
program" . 

The main program controls all the three sections of the program . 

The potential flow part uses six subroutines to determine the pressure 
distributions and the airfoil shape, including the separated streamlines. 
The viscous flow routine 'BLAYR' calls for three subroutines 'C0NV','AFSL' 
and 'LEASQ' to determine the boundary layer displacement thickness ^ , 
momentum thickness and the shape factor H. It also calculates the shape 
of the augumented airfoil by adding the displacement to the original 
airfoil shape. 

The separated flow region calculations are included in the main 
program. This part determines the rear stagnation pressure and the 
bubble mass flows in conjunction with 'Potential flow program' and 'BLAYR'. 

A. INPUTS 

The program input sequence is as follows: 

1. Main parameters of the program: The panel details and all the requi- 
red flow parameters. 
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2. The station distances for the ordinates of the airfoil. 

3. The airfoil thickness distribution; defined as Zj = (Zy-ZL) 

2 

where Zy and Zy are the upper and lower co-ordinates from the mid-section 
line at the given airfoil stations. 

4. The airfoil camber distributions; defined as Zq = ( Zy + Zy) 

2 " 

5. The panel widths for the whole range. 

6. The slopes of the upper surface of the deflected airfoil at the control 
points, calculated by a separate program from co-ordinates and 

7. The slopes of the lower surface of the deflected airfoil at the cont- 

< 

rol points , calculated as in 6. 

8. Specified pressures for points after the separation point if a separ- 
ation point is assumed to start the calculations. 

9. The angle of attack of the airfoil. 

B. OPERATION 

A diagram of the computer logic flow is shown in Fig. 20. 

The 'Potential flow program' prepares the airfoil for the solution by 
locating the panels and the control points. The source distribution and 
the vorticies are placed on the panels and at control points respectively. 
The solving of the simultaneous equations to determine the vortex strengths 
is performed by a standard IBM subroutine 'SIMEQ', which uses the Gauss 
elimination method. The output from the 'Potential flow program' is in 
the form of pressure coefficients at the control points and the new 
ordinates at the panel beginning points. The velocity distribution and 
the airfoil coordinates are used as inputs to 'BLAYR' for the determin- 
ation of the boundary layer characteristics and the augumented airfoil. 
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The calculations are done separately for the upper and lower surfaces. 
Since the ordinates of the airfoil are obtained at panel beginning points 
and the velocities are calculated at panel control points, subroutine 
'CONV converts the ordinates to the panel control points by a three 
point parabolic interpolation. 

The boundary layer displacement thickness determined by Head's 
entrainment method is added to the previous ordinates of the airfoil 
at the control points. 

The separation location is determined by comparing the values of H 
with a specified H5^p value. The control point at which H first exceeds 
H$£p is taken as the separation point at which separation will occur. 

The program is designed, however, to move the separation point only one 
panel at a time to avoid numerical instabilities. This is continued until 
the assumed separation point reaches the true separation point determined 
by the H distribution. The separation point is moved downstream if the 
H distribution fails to reach the specified value of H5£p. This allows 
free movement of the separation point depending on the pressure distri- 
bution. 

The separation point location determines the separation pressure, 
which is the value of the pressure coefficient at the separation' location 
chosen from the previous pressure distribution on the airfoil. The pres- 
sure on the trailing edge panel control point station is also set to the 
separation pressure to satisfy the Kutta condition. 

The boundary conditions are rearranged by specifying the separation 
pressure for the upper surface panels downstream of the separation panel 
up to the last panel on the airfoil (i .e. , trail ing edge panel) and a para- 
bolic increase from the separation pressure to the rear stagnation 
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pressure on panels up to the last panel of the flow field. The panels 
downstream of the trailing edge station on the lower surface of the 
airfoil are also specified by the same parabolic increase. The boun- 
dary conditions on the remaining attached flow panels are specified by 
the new slopes. The new slopes are calculated by the subroutine 'AFSL ' 
from the augumented airfoil surface points by matching a three point 
parabola. 

The potential flow routine now recalculates the new pressure 
distribution of the modified airfoil. This is continued iteratively 
until convergence is reached, i.e., the variation of the pressure is 
within 0.01. This is denoted by 'Convergence 1' (Fig. 21 ). 

After achieving 'Convergence 1' the program now calculates the 
separated region conditions. The location of the virtual origin of jet 
mixing flow (xQ,yo) and the mass flows entering and leaving the bubble 
are determined. The two mass flows and should be equal. Mass flow 
equality is acheived by iterating for the proper stagnating streamlines, 

R and R' , of the reattachment point (see Fig. 18). This is denoted by 
'Convergence 2'. The new rear stagnation pressure is used to recalculate 
the pressure distribution on the airfoil and the separated region. That 
is, we return to the potential flow and viscous routines. After 'Conve- 
rgence 1' is again reached the new stagnation pressure is found by 
satisfying 'Convergence 2'. This is iteratively continued until the- 
changes in the values of the reattachment pressure coefficient are with- 
in 0.001. This is designated by 'Convergence 3.' 

The last step in the program is the calculation of the recircula- 
ting mass flow in the separation bubble. The streamwise flow and the 
reverse flow masses are calculated and, if not equal, the value of Hgep 
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is changed so that the separation point is moved upstream to create the 
changes necessary for the mass balance (see Fig. 21) . 

This has to be iteratively continued until convergence is reached. 
This iteration involves all the other convergences (see Figs. 20 and 21). 
When this is reached the output corresponding to this iteration is the 
solution of the program. 

C OUTPUT 

The output from the computer are the following: , 

1 . The panel positions. 

2. Pressure distributions on both the upper and lower surfaces at the 
control points 

3. The ordinates of the airfoil and the separation bubble at the panel 
beginning points 

4. The position of the control points 

5. The slopes of the upper and lower surfaces 

6. Velocity distribution on the upper and lower surfaces at the control 
points 
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D. NUMERICAL INSTABILITIES AND THEIR REMEDIES IN THE OPERATION OF THE 
PROGRAM 

The boundary layer program exhibited numerous instabilities during 
its initial development. The cause of these instabilities could depend 
on a confluence of the various small irregularities in the different dis- 
tributions of the physical parameters. The values of H on the upper sur- 
face near the leading edge were sometimes exceeding the H^£p value and 
thus spread the instability to the downstream points. The same instabil- 
ity was also experienced on the lower surface near the stagnation point. 

An upper limit and lower limit for the H value were specified in order to 
damp out these oscillations. In addition to this, the comparison for the 
H value with the value was started only after about 15% chord to 
guard against the indications of , premature separation. The Felsch ex- 
pression for skin friction caused the program to find an exponent of a 
real negative number beyond a value of H = 3.0. Hence the upper limit 
was set at this value of H. 

During the development of the viscous flow program numerous inst- 
abilities were encountered. Any abnormally high values of s*, which caused 
wild variations were smoothened out by forbidding negative gradients of 6* 
for the first half of the airfoil. In spite of this, variations in 6* 
caused severe instcbil ities in the Cp distributions due to erroneous slopes. 
The 6* distributions were further smoothened by averaging the values with 
the values of previous viscous calculation since it was noticed that 
the variations were subsiding very slowly as the program went through 
the iterations. 
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Satisfying the Kutta condition by forcing the separation pressure 
on to the trailing edge point on the lower surface of the airfoil caused 
fluctuations of pressure distribution in the trailing edge region. This 
was confined to only the last 10 points. A second order least squares 
fit was used to smooth the calculated slopes in that region, which cured 
the instability. 

E. COMPUTER TIME AND COST ESTIMATES 

The computer runs were made on the IBM 360/44 and 370/145. The 
estimates are given in the following table. 


a 

Computer 

Time 

Total Number 
of Iterations 

Amount 

18.4° 

40 mins. 

72 

$ 120.00 

16.4° 

21 mins. 

44 

$ 68.00 

14.4° 

14 mins. 

27 

$ 36.00 


Initially a trial run with 69 smaller panels was attempted with 
the hope that the instabilities with 49 larger panels would be avoided. 
But the length of time required and the failure to cure the instabilities 
caused a return to the larger size panels. 

A typical convergence pattern and history for the case a = 18.4°, 

M = 0.135, RN = 2.2 x 10® is shown in Figs. 22, 23 and Table 2. 
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VI. RESULTS AND DISCUSSION 


The present analysis was used to find the pressure distributions 
on a GA(W)-1 17% airfoil at three angles of attack (18.4°, 16.4°, 14.4°) 
for which experimental data were available (Ref. 19). Fig. 25 shows the 
6A(W)-1 airfoil in the a = 18.4° position with the separation stream- 
line. The results of the pressure distributions are shown in Figs. 26,. 
27 and 28. It can be seen that the agreement with experimental data 
is good in the forward position of the airfoil and the separation pres- 
sure is predicted quite accurately. The position of the separation 
point is a little aft of the experimental value. The rear stagnation 
pressures also agree well with experiments. The values of the rear 
stagnation point pressure for separation and the separation pressures 
are given in Table 1 . 

The stagnation point pressures are close to the experimental data 
and H-separation values are also near those found in Ref. 2. 

The empiricism in this method is confined to (a) the position of 
the rear stagnation point, which is assumed from experimental data and 
appears to have slight influence on the results, (b) the well-established 
GOertler jet spreading rate parameter a, and (c) u^, the reverse flow 
velocity in the separation bubble. The assumption of a constant pres- 
sure separated region up to the trailing edge is such a well-established 
result from experiments that it can be considered to be a fact. 

The method does not restrict any of the separation variables. 

The separation point is determined by the boundary layer analysis for . 
each of the changed pressure conditions and is allowed to move. The 
mixed boundary condition program retains the separation pressure at 
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its specified value during the potential flow calculation, but once the 
separation point is altered, the sepration pressure changes to the value 
of the pressure at the new separation point. Thus there is no restric- 
tion applied on either the position of the separation point or its pres- 
sure. 

The empirical assumption of the position of the rear stagnation 
point is based on carefully measured experimental (Ref. 2 ) data,, and 
the results of the electric analog also showed that small movements of 
the rear free stagnation point did not produce any change in the stream- 
line pattern around the airfoil. 

Finally, the assumption of a uniform velocity u^ for the reverse 
flow was also based on carefully measured experimental data. The axi- 
symmetric separation model of Green (Ref. 15) also assumes a constant 
velocity reverse flow behind the base. 

The potential flow program calculates high positive pressure 

on the panels upstream of the front stagnation point on the lower 

surface. The panels up to the stagnation point were not included in 

★ 

the viscous analysis for this reason. The lower surface 6 was assumed 
to be zero for all the upstream panels since the flow is accelerating 
very fast in the forward direction. This avoided the problem of 
having pressure coefficients higher than 1 on the lower surface. 

The choice of u^,, though based on experimental data, turned out 
to be the key parameter for the final convergence of an acceptable 
solution, since u^ is used to find the reverse mass flow. Values of 
Ur/Ue = 0.12, 0.15, 0.18, and 0.20 were tried and the ideal value for 
which the solutions matched well with experimental data was Uj./Ug = 0.20 
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for both a = 18.4° and a = 14.4°. This indicates that the value of 
0.20 may be a universal term. 

The convergence for a = 16.4° was not achieved at u^ = 0.20 
but rather at a value of u^ = 0.22 u^. From the results of measurements 
on GA(W)-1 wings in Ref. 2 it was noticed that a = 16.4° represents the 
stalling condition (i.e.j the condition) and this is a highly un- 

stable condition. Measurements of the separated region data were diffi- 
cult to make due to the unsteadiness of the flow. We could perhaps 
attribute the present difficulty in the convergence at the same value 

of u to this unstable mode. 

r 

The MBC potential flow program used singularity distribution on 
the chord line to represent the airfoil. This would accurately predict 
results for thin airfoils, but for airfoils with large thicknesses, such 
as GA(W)-1, therecanbe significant errors in the pressure distribu- 
tions. A better potential flow solution which uses the singularity 
distributions on the surface could improve the theoretical solution to 
match more accurately the experimental pressure data. 

i 

Head's boundary layer method is an integral method and thus is 
limited by its assumptions. A more accurate boundary layer method 
could also improve the present solution. 

The accuracy of the method also depends upon the panel size. 

The pressure distributions shovy certain deviations from the experi- 
mental values at the points where there is a sudden change in the 
panel size. The choice of the panel size can greatly change this. 

By choosing a larger number of panel the accuracy of the method could 

be improved at the cost of increased computing time. 

Typical H and 6* distributions for the GA(w)-l airfoil are 

shown in figures 29 and 30. 
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Despite these limitations, this model has been shown to include all 
the significant physical features of separated wing flow. A computational 
method has been formulated which gives good surface pressure results. 

The one sensitive empirical value, u^/u^, may be a universal value of 

0.2; it is recommended that computation be performed for other wings to 
determine whether this is true. 

RECOMMENDATIONS FOR FURTHER WORK 

1. This program should be used for computing pressure distribu- 
tions on other airfoils for which separated-flow data are available. 

If data can be obtained for the separated region velocity field, the 
Ur/ue = 0.2 assumption can be tested directly; otherwise, the value 
will have to be inferred from the pressure field which results from 
various ur/ue values. 

2. An extension to the present model and program should be made 
to permit drag estimations. This could be done by computing pressure 
and velocity values on a vertical plane at the trailing edge station, 
and carried to a distance sufficient to insure that the flow is undis- 
turbed. The momentum deficit method could then give a drag value based 
on the detailed pressures and velocities at this plane. 

3. An improved mixed-boundary-condition potential flow program 
and boundary layer calculation method should be sought and mated to 
the separated flow model. 

4. Attempts should be made to extend this model and method to 
apply to multi-element airfoils, i.e., wings with flaps. 
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Figure 1 ga(w)-1 airfoil coordinates 


FIGURE 2 



STREAMLINE PATTERN-NAVI ER STOKES SOLUTION 

(from ref. 8 ) 
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Separation Point, Single Source or Source Distribution. 



Figure 3 simulation of separated region by sources 
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Comparison of Calculated and Experimental Data on 
NACA 63^-018 at a = 16° 


y/c 



Calculated Flow Pattern of Separated Flow, 
NACA ei^-OlS at a = 16°. 


Figure 4 source distribution model (from ref4 ) 
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ill U II 
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Figures equivalent airfoil model for separated flow (from ref. ii) 




Pressure distribution for flow with a separated 
wake (Jacob and Steinbach, 1974). 



Sketch of the new Jacob dead air flow region 
model with unsymmetric separation (Jacob 1975). 

Figure 6 Jacob's modified model (from ref.i.o) 
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Figure 7 experimental velocity plot, ga(w)-1 airfoil 
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Figure 12. details of the separated region 
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Figure 13 


EXPERIMENTAL SETUP FOR ELECTRICAL 
ANALOGY 
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Figure 14 streamline pattern from electrical analogy 
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Figure 18 geometry of the "inside problem 
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Figure 19 Matching of the existing flow with korst's flow at the separation point 
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Figure21 computer program convergence scheme 
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O STREAMWISE FLOW 



64 


Figure 22 typical plot showing bubble mass flow convergence (converge 
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Figure 23 separation Cp variation and convergence 
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Figure 2^ 


VELOCITY PROFILE IN THE SEPARATION BUBBLE 
ON THE UPPER SURFACE 
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Figure 26 pressure distributions, ga(w)-1 airfoil 
“= 18.4°, RN = 2.5 X 10‘, 





POTENTIAL FLOW 



Figure 28 pressure distributions, ga(w)-1 airfoil 
a = 14.4, RN= 2.9 X 10‘ 
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Figure 29 typical h-distribution on upper surface 
6AOD-1, a = 18.4°, RN = 2.5xl0^ 
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TABLE 2 


Typical Convergence History 



a = 18.4° 

, M = 0.135, RN 

= 2.2 X 106 



Iteration at Which Convergence is Achieved 


%EP 

Convergence 

1 

Convergence 

2 

Convergence 

3 

hsep 

91.75 




2.2 

86.75 

76.75 

3 



2.2 

71.75 

4 




66.75 

5 




71.75 

6 




71.75 

7 




66.75 

8 




61.75 

9 




61.75 

10 




61.75 

11 




61.75 

12 

12 



61.75 

13 

13 

13 


61.75 

14 

14 

14 

2.15 

56.75 

15 



2.10 

61.75 

16 




61.75 

17 




61.75 

18 

18 



61 .75 

19 

19 

19 



56.75 20 2.05 

56.75 21 

56.75 22 22 

56., 75 23 23 23 


51 .75 

55.75 

56.75 

56.75 

51.75 

56.75 
56.75 
56.75 
56.75 

24 

25 

26 

27 

28 

29 

30 

31 

32 

27 

31 

32 

32 

2.00 

51.75 

33 



1.95 

46.75 

34 




51.75 

35 




51.75 

36 




51.75 

37 

37 



51.75 

38 

38 

38 



1.90 
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TABLE 2 - continued 


46.75 

46.75 

46.75 

39 

40 

41 

40 

41 

41 

41.75 

42 



1 .85 

41 .75 

43 



41.75 

44 

44 


41.75 

45 

45 1 

EH 



5 



A 

ALFA 

ALFB 

B 

BETA 

C 

CF 

CPINV 

CPL 

CPR 

CPU 

DELSL 

DELSU 

DELX 

DHDX 

DQDX 

DZDX 

EL 

EN 

ETAM 

F 

FID 

FID! 

FIDL 


COMPUTER PROGRAM FORTRAN VARIABLES 


Array for coefficients of set of equations 
Angle of attack, a 
Tan a 

Array for coefficients of set of equations 

e = 1 - M2 

00 

Array for coefficients of set of equations 

Coefficient of friction 

Pressure coefficient of inviscid flow 

Pressure coefficient on lower surface of airfoil 

Pressure coefficient of rear, free, stagnation point 

Pressure coefficient on upper surface of airfoil 

Displacement thickness - lower surface 

Displacement thickness - upper surface 

Distance from the jet mixing "virtual origin" 

dH/dX 

de/dX 

dZ/dX 

Panel width 

Exponent of boundary layer profile 
n^p, non-dimensional ized coordinate in jet mixing 
Head's entrainment parameter in boundary layer 
4>d = Ud/Ue 

(j)d for the previous iteration 
‘J’d[_ “ '^d/^e lower surface 
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FMDOT 

Mp; streamwise component of recirculating mass 
separation bubble 

flow in the 

GAM 

Vortex strength 


HL 

Shape factor, H, on lower surface 


HSEP 

Specified H value for separation 


HU 

H values for upper surface 


IIUD 

Korst mass flow integral for the R streamline 


IlUJ 

Korst mass flow integral for S streamline 


IJ,IK, 

IL,IM, 

IN 

Counters 


IPRINT 

Printing counter 


ISEP 

Assumed separation location 


ISEPT 

True separation location 


ITERN 

Iteration counter 


lU 

Counter 


KOUNT, 
KOUNT 2 

Counters 


MDOTDL 

Mass flow into the separation bubble 


MDOTDU 

Mass flow out of the separation bubble 


N 

Total number of panels 


NA 

Number of panels on the airfoil 


NJ 

Panel number of the location of forward stagnation point on 
lower surface 

NL 

Number of panels on lower surface 


NT 

Number of airfoil points at which ordinates are 

given 

NU 

Number of panels on upper surface 


NX 

Total width of the arrays 


Q 

Height for reverse flow in separation bubble 
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R 

RMDOT 

SEPX 

SIGMA 

TCPR 

TEMPDL 

TEMPDU 

TEST 
TEST 1 
TEST 3 

TGAM 

THETAL 

THEATU 

TRIAL 1 
TRIAL 2 

UD 

UINF 

UL 

UNIT 

UU 

WL 

WLB 

WTE1 

WU 

WUB 

XA 

XC 


Recovery factor (not used for low speed flow) 

Reverse flow component of recirculating mass 
The x-coordinate of the separation point 
a, jet spreading rate parameter 
Temporary storage for CPR 
Temporary storage for DELSC 
Temporary storage for DELSU 

Convergence criteria parameters 

Vortex strengths (obtained from solving the simultaneous 

equations) 

Momentum thickness on lower surface 
Momentum thickness on upper surface 

Convergence criteria parameters 

Velocity of the streamline stagnating at the free stagnation 
point 

Free stream velocity 

Lower surface velocity 

Reynolds number 

Upper surface velocity 

Slope of the lower surface 

Input slope of lower surface 

Initial airfoil slope 

Slope of the upper surface 

Input slope of the upper surface 

X-coordinate of airfoil points 

X-coordinate of panel control points 
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XCPT 

XE 

XF 

XM 

XMI 

XO 

XR 

XTE 

YO 

ZC 

ZEL 

ZEU 

ZL 

ZT 

ZU 


X-distance of control point from panel beginning 
X-coordinate of a panel beginning points 
Final computation point 
X-coordinate of panel middle points 
Mach number 

X-coordinate of virtual origin (jet mixing) 

X-distance of rear stagnation point from trailing edge 
Trailing edge X-coordinate 

Z-coordinate of the virtual origin (jet mixing) 

Camber distribution 

Lower surface Z-ordi nates of airfoil at panel beginning points 

Upper surface Z-ordi nates of airfoil at panel beginning points 

Lower surface Z-ordi nates at airfoil points 

Thickness distribution 

Upper surface ordinates at airfoil points 
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COMPUTER PROGRAM INPUT AND OUTPUT DESCRIPTION 


A. The program input sequence is given below: 

(1) TITLE CARD 

The first card contains the title in columns 1-80. 

(2) PANEL DETAILS CARD 

This gives the total number of panels and its distribution 
in the field. 

N, NA, NT, NU, NL, IK (=1) are read in a (1015) Format. 

(3) FLOW PARAMETERS CARD 

This inputs the Mach number, Reynolds number. Recovery 
Factor (not used for low speeds) and some relevant positions. 

XMI, UNIT, R are read in a (3E20.6) Format and XTE, XO, XF, 
XCPT, WTEl are read in a (7F10.5) Format. 

(4) AIRFOIL COORDINATES CARDS 

The airfoil station positions, thickness distribution and the 
camber distribution are read successively. 

XA, ZT, ZC are read in a (7F10.5) Format. 

(5) PANEL LENGTH CARD 

The panel length on the airfoil as well as in the wake are 
given by this card. 

EL, is read in a (7F10.5) Format. 

(6) AIRFOIL INPUT SLOPES CARDS 

The airfoil slopes determined by an auxiliary program are pro- 
vided. The slopes correspond to the airfoil at the prescribed angle 
of attack. 

WUB and WLB are read in successively in a (7F10.5) Format. 
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(7) ANGLE OF ATTACK CARDS 


At present only one angle of attack can be read, but the program 

is designed so that it can be altered to read different angles of 

attack to get the results for various angles. 

IN (=0) is read in a (15) Format. 

ALFA is read in a (6F10.5) Format. 

Note: i) The value of ur/ue = 0.2 is included in the program. 

If convergence is not achieved this value may have to 
be changed. 

ii) The values of ZEU(2) and ZEL(2) have been provided. 

B. The program output is as follows: 

The first iteration gives the inviscid flow results. The viscous ef- 
fects are introduced in the second iteration. The value of NPRINT deter- 
mines the printing sequence. At present the first and second iterations 
are printed and every iteration corresponding to convergence 3 is printed. 
The iteration corresponding to convergence 4 is also printed. 

The details of the output are as follows: 

(1) Station points defining panels. 

(2) Ordinates of upper surface of the airfoil at the angle of attack. 

(3) Ordinates of the lower surface at the angle of attack. 

(4) Station points at panel midpoints. 

(5) Pressure coefficients on the upper surface at panel control 
points of the airfoil, and the wake. 

(6) Pressure coefficients on the lower surface at panel control 
points of the airfoil, arid the wake. 

(7) Slopes of the upper surface of the airfoil at angle of attack. 

(8) Slopes of the lower surface of the airfoil at angle of attack. 

(9) Velocity on upper surface of airfoil at panel control points. 

(10) Velocity on lower surface of airfoil at panel control points. 
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C SEPFLU/- - A PROGRAM FUR SEPARATED FLOW CALCULATIONS 

C FOP AIRFOILS AT LOW SPEEDS. 

C 

C THI:> PROGRAM COMPUTES THE SHAPE AND PRESSURE D I S TR 1 3UT I ON OF AN 

C. AIRFOIL WITH TRAILING EDGE STALL, INCLUDING THE SEPARATION tiUBBLE 

C 

C WRITTEN BY SHAKAD N . NA IK 

C WICHITA STATE UNIVERSITY 

C APRIL 8, 1977 

C 

C THE POTENTIAL FLCW PART CONSISTING OF SUBROUTINES 

C FLSOLV, SETUP, SLOPE, PARINT AND PARINB WAS PROVIDED 

C BY MCDONNELL-DOUGL AS AIRCRAFT COMPANY. 

C 

C THE SEQUENCE IN ^HICH THE DATA WILL BE READ IS PROVIDED 

C IN REPORT AR77-2, DEPARTMENT OF AERONAUTICAL ENGINEERING, 

C WICHITA STATE UNIVERSITY, WICHITA, KANSAS 6720R 

C 

CJMMCN /READ/ A T I T L E ( 2 0 ) , N , N A , NT , X M I , UN I T , R , XTE , XO , XE , XCP T , W T E 1 , 

1 XA ( 51 ) , ZT (B 1 ) ,ZC( SI ) , EL 166) , ALFA, WUEND 

COMMON / I iN T/ I U , IL , I S , I LAM , I SE P, NJ ,,NL, NP 1 , NM 1 , N2 , NUl ,NL I , I J , I K 
COMMON /GEOM/ X E ( 66 ) , XM ( 6o ) , XC (66 ) , ZUB I 5 1 ) , ZL B ( 5 I ) , Z A { 66 ) , Z 3 ( 66 ) , 

1 ZEU (66) , ZEL ( 66) , ZZ( o6 ) , ZU(51 ) ,ZL( 51 ) 

COiMMCi'J /pan/ WUB(50J ,WLB(50) J^Cd(5D) ,WU(6t>) ,WL(66) ,BETA,ALFB 
COMMON /PAN2/ OP L ( 66 ) , CPL ( 66 ) , UU (66 ) ,UL I 66 ) , I TERN 
COMMON /MAT/ A ( 66 , 6 7 ) , B ( 66 , 66 ) , C ( 1 32 , 1 32 ) 

COMMON /MATR/ G A M ( 1 3 3 ) , S E ( 66 ) , TG AM { 1 33 ) 

COMMON /bLR/ HU( 50) ,HL(50) ,THcTAO( 50), THETAL(50) ,0ELSU(50), 

1 UELSL(50) , S E!^X, HS EP., CP I NV{ 50) ,NJ 

C 

REAL MCGTOO, MDOTDL , I 10 J , I 1 UD , Y BU ( 6 6 ) 

READ IN BASIC AIRFOIL DATA AND SET UP PANELS 

READ (5,150) (ATITLEd ) ,1 = 1,20) 

READ (5,155) N,NA,NT,NU,NL,IK 
READ (5,160) XM I ,UN I T, R , XTE , XO ,XF, XCPT , WTEl 
C 

C?R = 0. 1 
HSEP=2.2 
SEPX=1 00. 

IS£P=NA 
K JUN T2 = l 
SIGMA= 12 
IPRINT=0 
IIUJ=0.39Q 
ETAM=0 . 399 
UINF =0. 135--;d 085. lOA 
SUH=C. 0 
N01=KU+1 
NL 1=NL+1 
N A 1 = N A + 1 
I TER.\= 1 
C 

DiJ 5 I = 1 , N 
OPU( [)=o.j 
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CPU n=o.u 
5 CUNFINUE 

READ [i\ AIRFCIL Cu'GRDINATES AWC SLOPES 

READ {5,lfiS) (XA(I)tI = l»M) 

READ { 3,165) ( ZT(I ) ,I = l,.\T) 

READ ( 5, 165 ) ( ZC N ) , 1= 1, ^T) 

READ (5,165) ( E L ( 1 ) , I = i , N ) 

READ (5,165) ( WU 3 ( I ) , I = 1 , iMU ) 

READ (5,165) ( WL b ( I ) , I = 1 , NL ) 

READ ( 5,165) (ZUd ) ,I = 1,M) 

READ (5,165) ( Z L ( I ) , I = 1 , M ) 

READ ( 5 , 17C, ENU= 145) I M, ALFA 

ALF6=ALFA/57 .29578 
ALFb = rA,'J ( ALFb) 

XE ( 1 ) = 0 .0 

DU 10 U1,N 

XE( I + l ) =XE ( I )4-EL ( I ) 

XM( I ) = XE( I ) 4-0. 5 -EL ( I) 

XC ( I )=XE( I ) 4-XCPT*EL ( I ) 

10 COMTIOUE 

KUUN r=i 
CALL SETUP 

00 15 I=1,N 
WU( I )=0.0 
aL ( I )=0.0 
UU( I )= J.O 
UL ( I 1^0 ,0 
GAM( I )=0.0 

15 CUNT i:\lUE 

DU 20 I=1,NP1 
ZA( I ) = 0.0 
ZEU ( I ) =0 .0 
ZEL ( I ) =C, 0 
GAM(,N4- I ) = 0 .0 
5E ( I )=0.0 
20 CONTINUE 

ISEP=0 
l'J = 0 
IL = 0 
I;'1=0 

1 J=0 

CCHEK = 0 .0 


rmiALIZE THE SLOPE VALUES 


CALL SLOPE 
ZA( 1 ) = ZUb( 1 ) 
Z6( 1 )=ZL5( 1) 
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OJ 2 5 I = l,:vJA 

ZA( I + l ) = Z'Jd( I + l )-XE( I-H ) «ALFd 
Z3( H-l)=ZLtt( I+L)-aE(I + 1J *ALEB 
25 CGAiTINUE 

SE( 1 ) = 2.0^(WTE1-ALFB) 

N=i\A 
ISEP=NA 

IF ( ITEiXN. EG .1 ) GO TO 40 
.30 C?L( NL )=CPU( ISEP) 

Du 3 5 I =iN A 1 f N 

CPU! I ) =CPU( I SEP) +(oPR-CPU( ISEP) J (3.*(XC( I )-100.) /( lOO.-SEPX) )=s 
1 * 2 ) 

CPU I ) =CPL ( NL) + ( CPR-CPL ( NL ) ) ^--I (3 XC ( I ) -1 CO . ) / ( 100 .-SEP X ) ) ««2 ) 
35 CONTI.NUE 

CCOPLTE AIRFCIL PRESSURE DISTRIBUTION AND BUBBLE SHAPE 
40 CALL FLSCLV 

IF { I TER.'J.GT .1 .AND. I L. EQ. 1 ) GU TO 50 

PRINT TITLE AND RESULTS FCR THE PRESENT ITERATION 

45 WRITE (6,175) ( A T I TL E ( I ) , I = 1 , 2 0 ) , A L F A , I T ERN 
WRITE (6,130) ( XE( I ) ,1 = 1 ,NP1) 

WRITE (6,135) (ZEU( I ), I = l,NPl) 

WRITE (5,193) ( ZEL( I ) , I = 1,NP1) 

WRITE (6,205) ( XM ( I ) , I = 1 ,N ) 

WRITE (5,190) (CFU( I J , I = 1,N) 

WHITE (6,200) (CFL( I ) , 1 = 1 ,N) 

>(RITE (6,210) (XC( I ) , I =1 ,N) 

WRITE (6,2 15) ( WU( I ) ,1=1 ,N) 

WRITE ( 6,220) (WL( I ) ,1 = 1, N) 

WRITE (6,225) (UU( I ) , 1 = 1 ,N) 

WRITE (6,230) ( U L ( I ) , I = 1 , N ) 

C 

IF(IPR1M. EG. 2)ST OP 
50 IF ( IPR I WT.EQ. I ) GC TO 135 
C 

IF ( I tern. EG. 1) GO TC 55 
GO TC 70 
C 

5 5 00 6 0 I = l,i\A 

CPINV( I )=CPJ ( I ) 

ZU ( I ) =ZEU( I ) 

ZL( I )=ZEL( I ) 

6 0 C U N T 1 N U E 
r 

NAl=NA+l 

WRITE (6,130) (XE( I ) ,I = 1,NA1 ) 

WRITE (6,135) ( ZeU( I ) , 1= 1,NA1) 

WRITE (6,195) ( ZEE ( I ) , 1 = 1 ,NA I ) 

C 

ISEP=NA-4 
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SEPX^XC ( I SEP ) 

NU=ISEP 

XR = ( lOJ.-ScPX) /3 . 
iMR = XR+ L 
I'J R +i\ A 

I reRN=I FERN+l 
ISEP1=ISEP+1 
C 

IJ J 6 5 I = I S E P , N A 
65 CPU( I ) = CPU( ISEP) 

GG TC 3C 
C 

CALCULATE THE BOUNDARY LAYER CHARACTERISTICS 
70 CALL 3LAYR 

IE ( I rERN.GT,3 .AND. IL.EQ.l) GO TO 75 


PRINT BOUNDARY LAYER CHAR AC T ER I S T I C S FUR THE CURRENT ITERATION 


WRITE 
WRITE 
irtR I T t: 
WR I I E 
WHITE 
WRITE 


{ 6,235 ) ( HOI I ) ,1 = 1 , I SEP) 

( 6,240 ) ( HL( I ) , 1=1 ,NL ) 

(6,2451 ( THE TAU( I) , 1 = 1 , I SEP) 

(6,250) ( THETALI I ) ,I=L,NL) 
(6,255) ( JELSU( I ) , I = 1,ISEP) 

(6,260) (DELSL ( I ) , I =I ,NL ) 


75 IL=L 

XR= ( 100 .-SEPX) /3 . 

NR=XR+ 1 
N=NR+NA 

IF (\.Gr,63) STOP 
IF (KUUNT.EQ.l) GG TO 00 
GO TO 8,5 


8 0 TK IAL1 = CPJ( I SEP- 3) 
KGUiNT = KGUi'iTHl 
GO TC 3C 

S5 TRIAL2=CPU( ISEP-3) 


PRINT TRIAL PRESSURE COEFFICIENTS 

WRITE (6, 88) TR lALl, TRIAL2,CPU( ISEP) 
38 FORMAT! 50X, 3F 10 .5 ) 

TEST = TRIAL 1- TR I AL2 

IF { ABS ( TEST ) .LT .0 .0 1 ) GO TO 90 

TR I AL1 = CPU( I SEP-3) 

K0UNT=I<GUNT^-1 
GO T i_i 3 0 


90 ISED1=ISEP+1 


SEPARAflON BUBBLE "MASS-IN / lASS-OJT" CONDITION 


CU( I SE^ ) =UU( ISEP )"‘C I i\F 
UL (NA)=UL( NA)-UINF 
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XO=iO.*THE TAU( I SEP ) 

YC=CELSU( I SEP) H PET AU( 1 SEP) 

C 

OU 9S I=ISEP,N 

9S YbU( i )=0ELSJ( ISEP)-YU+{ (XC( I )- SE PX ) «EF AM/S I GMA ) 

C 

E.'j=2./(HL (NL )-l) 

UO=U INF^SQRT (CPR-CPU ( [SEP) ) 

F IO=UO/L:G( ISEP ) 

FIDl=0.fcL2 
TEST 1=0.0 
C 

100 IlUD = 0. ll402-0.2 0'-f5*FI0<-l.0317=!=FI0«FID 

MDOTOU= (UU( I SEP) *{ XE{ NAl )-XC ( I SEP) + XC) / S IGMA)=i= ( I LUJ- I lUO) 
MOUTOL=UL( NA) >:=DELSL( iMA) { 00/ UL ( NA ))«=;={ Ei\l+ 1 .)) «EN* (( EN + 1 .)«*(- 1 . 
1 /LN) )/(EK+l. ) 

PRINT THE MASS FLOW IN ANO MASS FLOW OUT 

rtRITEIo , 102) MDL'TOU.MDOTCL »UD 
102 FORMA n 1 X, jF 15 .5 ) 

TEST=MOCTCU-MOaTOL 
IF lAdS( TEST).LE.O.Jai) GO TC 105 
HOLC=FIO 

FIO = FIO-TEST'-M F lOl-F I0)/(TEST1-TEST) 

FIC1=HULD 
TEST1=TEST 
UO = F IO*UU( I SEP ) 

GO TO 100 

PRINT REAR ATTACHMENT POINT PRESSURE 

10 5 CPR=CPU ( I SEP ) +• ( (U)/UNFJ««2) 

WRITE (6t265) CPR 
IF I KCON rZ .EG. 1 ) GG TO 110 
GU TO 115 

110 TCPR=CPR 

KGUNT2 = !<CUNT2-t- 1 
(; u T u 3 0 

115 TEST = aBS ( TCPR-CPR ) 

IF { TE ST.LE.O. 001 ) GO TC 120 
TCPR=CPR 
GG TC 30 

dJtibLE VORTEX HASS FLOW CG:mCITIUN 

120 ZUT=-C.C7 

FIOL=UD/UL (NA) 

OELX=XE { KAl ) +XQ-SEPX 

G=ZEU( NA ) 1-XC {.JA ) ^;=ALF6-Z0T<-Y')-( 2. 3 0 39 3-DELX/S IGMA )-DELSU( I SEP ) 

1 — JEL5L(NL)^(EN+i)^IF10L*^EN) 

FMUOT= ( UU( ISEP ) *'JELX/5 IGMA ) * ( I lUO-O . 23 2 1 +• 0 . 2 50 5 A ) 

R r-1 D G T = 0 - 0 . 2 0 U U ( I S E P ) <■ 0 . 0 3 8 6 * U 0 ( I S e P ) - 0 E L X / S I G M A 
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TEST3=AdS ( RMDUT )-ABS(FHUGT ) 


PRINT dJdbLE i'lASS FLU..^ VALUES 


/.RITE (6,122) RML:Gr,F,.MDCI,MSEP,SEPX,CPU( ISEP) 
12 2 FORMAT ( IX, 5F 15 .5) 


C 


C 

C 

C 

c 

c 

c 


r 

u 


c 

c 


IL= IL+1 

IF ( TESTJ.GT.0.0) GO TO 125 
GO TC 13C 


125 HSEP=HSEP-U.05 

UU(ISEP)=UU( I6EP)/UINF 
UL (MA ) =UL ( NA )/ J INF 
GO TC 30 

PRINT separation POINT SHAPE FACTOR 

125 WRITE(6,126) HSEP 

12 6 FORMAU I X, • hSEP = *,F5.2) 

130 IPRINT=1 
GO TC 45 
135 ISEP=IS£P+i 
IPRI i\T=2 

CPU( ISEP)=2^CPU( ISEP-l )-CPU( ISEP-2 ) 


00 140 I=ISEP,NA 
CPU! I ) = CPU( ISEP) 

140 continue 

GO TO 30 
145 STOP 


150 

FORMAT 

155 

FORMAT 

16 0 

FORMAT 

165 

FORMAT 

170 

FORMAT 

175 

FORM A r 

1 

2 

180 

FORMAT 

18 5 

FORMAT 

190 

FORMAT 

19 5 

FORMAT 

20 0 

FORMAT 

205 

FORMAT 

210 

FORMAT 

21 5 

FORMA r 

22 0 

FOR/! AT 

22 5 

FORMAT 

230 

FOR-IA I 

235 

FORM AT 

240 

FORMAT 

2^5 

format 


(20A4) 

( 1 i J I 5 ) 

( 3E20 ,6/ 7F10. 5) 

( 7F10.5) 

( I5,6F10.5) 

( IHl ,25X , ‘A IRFGIL nITH VISCOUS EFFECTS INCLUDING THE ', 
'PUSSIrflLITY UF TRAILING EDGE ST AL L ' /26 X , 20 A4/ 1 HO , 

•ALPHA =',F1J.5,* DEGREES' ,5X, • I TERATI CN NUMbER =',13/1H0 
{'JSTAriGN POINTS ■ OEF IN ING PANELS'/(1H ,10F12.5)J 
( 'OORDINATES OF UPPER S URF AC E • / ( IH ,10F12.5)) 

COPRESSURE COEFFICIENTS CN UPPER SURFACE'/(TH ,10F12.5)) 

I 'OORDINATES OF LCoER SURFACE'/(iH ,10F12.5)) 

COARESSUFE CJEFF IL IENTS on lower SURFACE'/! IH ,1CF12,5)) 
COSTAnOM POIMTi AT PANEL M I DP C I NT S ' / ( 1 H ,10F12.5)) 
CJSTATIGn points at panel control PCINTS'/IIH ,1UF12.5)> 
CUSLCPES CF UPPER SJRFACE'/(1F ,10F12.5)) 

('OSLGPES OF LOWER SURFACE* /(IH ,10F12.5)) 

CDVELuCITY ON UPPER SURFACE XC'/(1H ,10F12.5)) 

CUVELJCITY ON LOrttR SURFACE XC * / ( 1 H ,10F12.5)) 

(• hU', 10 FI 3. 5) 

{' HL',10F1J.5) 

( ' ThETAU ' , 1 JF 10. 5) 
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250 FORMAT (’ TH ET A L * , 1 JF I 0 . 5 ) 
255 FORMAT (' CELSU SIOFIO.S) 
260 FORMAT (• CELSL •,UFlu.5) 
265 FORMAT 15X,'».?R = ',F10.5) 
EMD 
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^UriRiJUr INt SETUP 


C THIS SLBRUUriTE PcAUS DATA DEFINING AIRFOIL SHAPE, SETS UP PANELS 

ALONG THE AIRFOIL CENTERLINE AND FINDS COMPUTING RATIOS OF PANELS 
THIS HAS TO BE DONE ONLY ONCE FOR EACH AIRFOIL. 

COMMON /READ/ ATITLEI20) ,N,NA,i\T,XMI,U!MIT,R,XTE,XO,XF,XCPT,WTEl, 

I XA(5I),ZT{51),ZCli>l),EL(66),ALFA,WUEND 

COMMON /IN T/ I U , IL , I S , ILAM, ISE P, NO ,NL, NPi, NMl, N2,NU1,NL 1 , IJ , IK 
COMMON /GEOM/ X E ( 66 ) , XM ( 66 ) , XC I 65 ) , Z UB { 5 L ) , Z LB ( 5 I ) , Z A ( 66 I , Z B ( 66 ) , 
I ZEU{o6l,ZEL{66),ZZ(66),ZU(5I),ZL{51) 

COMMCN /pan/ /^U5(3G) , ',^Li3 ( 5 0 ) , WCd ( S 0 ) , w U ( 66 ) , WL ( 6 6 ) , d ET A , AL F B 
COMMON /PAN2/ CPU( 66) , CPLI66 ) , U'J(66 ) ,UL{ 66 ) , ITERN 
COMMCN /MAT/ A ( 66 , 6 / t , B ( 66 ,66 ) , C ( 1 3 2 , 1 32 ) 


INDICATOR IK = 0 VISCCCS SOLUTIONS AT VARIOUS ALPHA 
IK = L INVISCID SOLUTIONS AT VARIOUS ALPHA 
IK = 2 VARIOUS INVISCID SOLUTIONS 

IF NU OR NL IS LESS THAN NA, SOLUTION WILL bE FOUND FOR THE MIXED 
BOUNDARY CONDITIONS OF GEOMETRY SPECIFIED UP TO PANEL NU OR NL 
AND PRESSURE COEFFICIENTS ON SUBSEQUENT PANELS. 

5 aRIIE (6,30) ( ATI TL6( I) ,1 = 1,20) 

wR I T E { Ci , 3 5 ) N , N A , N T 

WRITE (6,'+J) XMI ,UNI rtR,XTE,XO,XF,XCPr,wT£l 
WRITE (6,45) ( XA{ I) , I=1,NT) 

WRITE (6,50) ( ZTU ) ,I = 1,NT) 

WRITE ( 6,55) ( ZC(I ) ,l=l,NT) 

WRITE (6,60) (Wue(I),I=l,NU) 

/<RITE (6,65) ( WLE( I ) ,1 = 1 ,NL) 

PI=3 .1415927 
TPI=2.G*PI 
XM2 = XMI>!‘XM I 
BETA=SQRI( 1.G-XM2) 

NP1=N+1 

NM1=N-1 

N2=2.0«N 

FIND COMPUTING RAT US OF PANELS 

Du 10 I = 1 , N 
A{ [ ,NP1)=0.0 
DO iC J=1,N 
A( I , J ) =0.0 

B(I,J)=(L.O/TPI)*(ALOG(AuS{(XC(I)-XE(J+L))/(XC(I)-XE(J))))) 

10 CONTINUE 

OU 15 I=1,N 
DU 15 J=L,N 

TL = AlOG( A6S( ( XM( I )-XE( J+1 ) ) / (XM ( I )-XE( J) ) ) ) 

A( I , J) = ( r 1-1 .0 )/PI-Tl*(XM( I )-XE( J) ) / ( P I =!'EL { J ) ) + A ( I , J) 

A( I , J+-1 ) = l .0 /P I + T1* (XM ( I )-XE { J ) ) / ( P I’'--EL { J ) ) + A( I , J + 1 ) 

15 CONTINUE 
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CU 2 0 [ = 1 , M 
DO 20 J=1,NP1 
A( I , J ) =A ( I , J ) /oE TA 
20 CONTINUE 
C 

ZUB( 1 ) = 2T( 1 ) 

ZL3( 1 ) =ZT( 1 ) 

CALCULATE SLOPE AT THE CONTROL POINTS 
INTEGRATE THE SLOPES TO OBTAIN AIRFOIL ORDINATES 

CO 25 I=1,NA 
XX=XC( I ) 

CALL PAR.INT ( X A , ZC , 1 , NT , XX , X X , Z , 1 , Z I ) 

WCB{I)=W1 
XX = XE( I + i) 

CALL PARIN3 (X A , ZT , 1 , NT , XX , X X , Z , W , Z I , I ) 

IF (I.EQ.l) Z=ZT(2 )--SQRT (XE(2)/XA(2) ) 

CALL PARINT ( X A , ZC , 1 , NT , XX , X X , Z L , W , Z I ) 

ZUB( I+l J=Z+Z1 
ZLBI I+1)=-Z+Zi 
25 Continue 

RE TURN 

30 FORMAT ( iHi,25X, • A IRFCIL /»ITH VISCOUS EFFECTS INCLUDING THE 
1 'POSSIBILITY OF TRAILING EDGE S TALL * / 26/ , 20AA/ LHO ) 

35 FORMAT ( I HO , 1 7 X , • TO TA L NCVBER OF PANELS ='»IA/ 

1 13X,»NUYbER OF PANELS IN AIRFOIL -',14/ 

2 7X, 'NUMBER OF POINTS DEFINING AIRFOIL =',I4) 

40 FORMAT (29X,«MACH NUMBER =',F12.5/ 

1 24X,'JNIT REYNOLDS NO =',E12.5/ 

2 25X, 'RECOVERY FACTOR =',E12,5/ 

3 19X, ' TRA IL ING EDGE STATION =',E12.5/ 

4 hX, 'DISTANCE FROM T.E. TO BUBBLE CLOSURE =',E12.5/ 

5 L5X, 'FINAL COMPUTATION STATION =',E12.5/ 

6 9X, 'DISTANCE TO PANEL CONTROL POINT ^*,E12.5/ 

7 L3X, 'INITIAL AIRFOIL SLOPE = ' , E 1 2 . 5/ IHO ) 

45 FORMAT ('OSTATlOiN POINTS DEFINING A I R F 0 1 L ' / ( IHO , 10 F I 2 . 5 I ) 

50 FORMAT ('JuRDlNATES OF THICKNESS 0 ISTR I BUT ION ' / ( IHO, lOF i 2. 5 ) ) 
55 FORMAT ('OOROINAT6S OF CAMBER D I STR I bUT I ON ' / { 1 HO , 1 OF 12 . 5 )) 

60 FORMAT ('OSPECIFIED UPPER SURFACE SL OP 6 S '/( I HO , I GF 1 2 . 5 ) ) 

65 FORMAT I'jSPECIFIED LOViER SURFACE SLOPE S ' / ( 1 HO , I OF 1 2 . 5 ) ) 

END 
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SUBROUn^F SLOPE 
INITIALIZE THE SLOPE VALUES 

COiMMCN /READ/ A T I T L £ ( 2 0 ) » ,\ » N A , NT , X M I , UN I I , R , XT E , XO , X F , XC P T , W TE L , 

1 XA(S1),ZT(5L),ZC{51),EL{66) ,ALFA ,>UEND 

COMMON /INT/ lU, IL,IS, RAM, I S E P , NJ , NL , NP I , NM 1 , N 2 , NL I , NL 1 ,IJ ,IK 
CGMMGri /GEOM/ XE ( 6 6 ) , X M { 66 ) , XC ( 6 6) , ZUB (31 ) ,ZLB(5L),ZA(66),ZB(66), 
1 ZEU{66),ZEL{66),ZZ(56),ZU(51),ZL(51) 

COMMON /PAN/ WUBI50) ,WLB(30) ,WCri(50) ,WU(66) ,WL(66) ,8ETA, ALF3 
COMMON /PAN2/ CP U ( 66 ) , CP L ( 66 ) , UU ( 66 ) , UL ( 56 ) , I TE R N 
COMMON /MAT/ A(66,67) ,B(6o,66) ,C( 132, 132) 

ZEU ( 1 ) =ZT ( 1 ) 

DO 5 I=2,NU1 
ZEU{I)=ZOB(I) 

3 CGNT INJE 

ZEL ( 1 ) =ZT( 1 ) 

DO 10 1=2, NLl 
ZEL( I) = ZL3( I ) 

10 CONTINUE 

DO 15 I=1,NU 
rtUn ) = aUB( I ) 

15 CONTINUE 

00 20 t=l,NL 
rtU I )=wLB( I ) 

20 CONTINUE 

RETURN 

END 
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SU3RUUTI']E FLSOLV 


THIS oljBRLiUr liME TAKES GECMET«Y, ESTABLISHES PRESSURES UN SPEUIFIE 
PANELS and THEN SULVES FCR PRESSURES OR GEOMETRY AS APPROPRIATE. 
THE SHAPE AND PRESSURE DISTRIBUTIUN IS PRINTED OUT. 

CCi-lMLN /READ/ AT ITLE(20) , N , N A , NT , XM I ,UN I T , R , XTE , XO , XF , XCP T , WTE 1 f 
1 XAI 51) ,ZT(5L) ,ZC( 51) ,EL(66) ,ALFA,WUEND 

COMMON / INT/ lU , IL , IS , ILA.M, I SEP,NU ,NL, NP1,NM1,N2 tNUl ,NL 1 ,I J , I K 
CG.MMCN /GEOM/ XE (66) , X M ( 66 ) , XC ( 66 ) ,ZUB(51) ,ZLB{5l),ZA{66),ZB(66), 
1 ZEUI66) ,ZEL(66) ,ZZ(66) ,ZU(51) ,ZL(51) 

C'J.'^MCn /pan/ wU 6(50) ,V(Lri(50) ,WC6( jO) ,wU(66) ,WL(66) ,BETA,ALFB 
COMMON /PAN2/ CP L I 66 ) , C P L ( 66 ) ♦ UU ( 6 6 ) , UL ( 66 ) , I T ERN 
CUMMO.N /mat/ A( 66,67 ), 3( 66,66) ,C I 132, 132) 
common /MATR/ gam I 133) ,SE (66 ) , TGAM ( 133 ) 

DIMENSION wJI(66), ImLI(66), XCI(66), CPT(66) 

DIMENSION WCI66) 

NX=132 
DO 5 I=1,N 
GAM( I )=0 .0 
GAM( N+ I ) =0.0 
SE(I)=0.0 
5 CON T I NOE 

N2P l=2*N + l 
GAM( N2P1 )=0.0 
SEINPl )=0.0 

ESTABLISH BOUNDARY CONDITIONS ON APPROPRIATE PANELS 

T2=C.5*XCPT 
T3=0 .5«( 1 . 0-XCP T ) 

N2 = 2.0’SN 
NP1 = N+ 1 

\ 

DO 10 I=1,N2 
DO 10 J=1,N2 
C( I , J ) =0 .0 
10 CONTINUE 

DO 45 1=1, N 

IF ( ITERN.EQ.l ) GO TO 15 
GO TO 20 

15 IF ( I .G I. NO) GO TiJ 35 
GO TO 25 

2 0 IF ( I .GE..^JU) GO TO 3 5 
2 5 GAMI I ) = aU( I ) 

IF II.EO.l) G.aM ( 1 ) =GaH ( 1 )-T3*SE( 1) 

CO 3 0 J = 1,.N 
C( I , J) ^d( I , J ) 

30 CCNT INuE 

C( I ,N+ I )=T2 
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c 


c 


c 


c 

c 


c 


c 


c 

c 

c 

c 


c 


c 


u 


IF ( I ) C( I ,N + I-1)=T3 

GO TC 45 

3 5 CPO( n = 2.*n .J-S CRT ( I .0-CPiJ{ I) ) ) 

GAM ( I ) =CPU( I )-A ( I , 1 ) ( 1) 

Cl I , I )=-1.0/BETA 

00 4 0 JL=2,<NiPl 

J=i\ + Jl-1 

C( I , J) =A ( I , J1 ) 

40 CGiMT I .'i Jc 

4 5 CGi\TI\OE 



DO 

8 0 

1 = 1 

, N 




IF 

{ ITERN 

.EO.l ) 

GO 

TO 50 


GO 

TO 

55 




50 

IF 

1 I 

• G T • 

NL) GO 

TO 

70 


GO 

TO 

6 0 




5 5 

IF 

1 I 

• G E • 

i'JLl GO 

TO 

70 

60 

GAMl N+ I ) = 

'aLI I ) 




IF 

1 I 

.FQ. 

1 J GAM(N+1J=GAM(N+1) + T3*S£( IJ 


DO 

65 

J = 1 

,N 



65 

Cl N + I 

, JI = 

ol I , J ) 




C I i\ + I » N + I ) T 2 

IF (I.,ME.l) C(:N+I,N4I-1)=-T3 
GO TO ao 

70 CPL ( n=2 .^( I ,0-SGPT ( 1 . 0-CPL ( I ) ) ) 

GAM( N-t-I )=CPL ( I )-A( I , I ) *SE( L) 

C(N+I , I )=1 .O/dEFA 

DO 75 J1=2,MPL 
J=N+ Jl-L 

C(iM+I , J)=A( I ,J1 ) 

75 CONTINUE 
80 CONTINUE 

SOLVE FOR SOURCE AND VORTEX STRENGTHS 

CALL SIMGQ ( C , GAM , TGAN' , N 2 , NX ) 

DO 35 1=2, NPl 
SE ( I > = TGAM(iM + I-l ) 

85 CONTINUE 

DO 1 L 5 I =1 ,N 
I } = 0*o 

CPT ( I ) = A(' I ,cvi?l ) =!=SE (NP 1 ) . 

DO 90 J=i,N 

CPT( I )=CPT( I )*-Al I ,J)4'SE( J) 

9 0 RC ( I ) = NC ( I ) H3 ( I , J ) TGAM ( J ) 


U = I oE ( I ) + XCPT=!={ SE ( I + 1 )-S£ ( I ) ) ) / 2. 0 
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AU( I )=WC( I)+T1 
WL { I ) =wC ( I )-Ti 

UU( I ) = i .0-0.5«(CPT( I )-TC.Ai>) ( I )/ tJETA ) 

IF ( I TERN. NE . L ) GU TG 95 

IF (I.LE.NU) UU( n=lJUl I)/SQRT{ l.OfWUd )=!'*2) 

GO TG 100 
C 

95 IF (I.LT.iNfJ) UU{ I)=UU( n/SQRT( 1.0+WU( I )«-2) 

100 CP'J( I ) = 1.0-UU( I)’!'-P5/ABS<UU{I)) 

UL ( I ) = 1 .0-0. 5=!'( CPT I I ) + TGAi'^( I )/BETA ) 

IF ( ITERN.NE .1 ) GO TG 105 

IF (I.LE.NL) UL ( I ) =UL( I) /SORT! 1.0 + WL ( I ) «*2 ) 

GC TG 110 

105 IF (I.LT.NL) UL ( I )=UL( I ) /'SGRTI 1.0 + WL { I )**2 ) 

110 CPL(n = 1.0-UL(I)=!'«3/AriS(UL{I)) 

115 CONTINUE 
C 

G The AIRFGIL slopes can NOrt 8E INTEGRATED TO FIND 

C THE SHAPE CF THE AIRFGIL AND SEPARATION BUBBLE. 

C 

ZEU( 1) = ZT( 1) 

WUI I 1) =SE( 1) /2.0 
ZEL(1)=ZT(1J 
>--LI ( 1 ) =-SE(l )/2 .0 
XCI ( 1) =0.0 
C 

GO 120 1=2, NPl 
wui ( n = rtU( i-i) 
ali ( n =aL( i-i ) 

120 XCI ( I ) = XC( I-l) 

c 

ZEU(2)=1.325 
ZEl ( 2) =-2.25 
C 

GG 125 I =3, NPl 
XX = XE{ I ) 

XO=XE( I-l) 

CALL PAR I NT ( X C I , WU I , 1 , N P 1 , XO , XX , Z , W , Z I ) 

ZEU ( I ) =ZEU ( I-l ) + ZI 

CALL PARINT ( XC I , aL I , 1 , N FI , XO , XX , Z , A , Z I ) 

125 ZEL( I) = ZEL( I-1)+Zl 
RE TURN 
End 
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SUtiRCUTINE 8 La\YR 


THIS ROUTINE UUKiPUfES THE tiCUNDARY LAYER CHARACTERISTICS 
USING THE VELCCITIES FRO) FLSULV. 

COMMON /READ/ A T I T L E ( 2 G ) , N , N A , N T , X M I , U N I T , R , X T E , X 0 , X F , X C P T , N T E 1 , 

1 XA(51),ZT(dI),ZC(51),EL(66) ,ALFA , aUEND 

COMMON /INT/ lU , IL , I S , ILAM, I SEP ,NU ,NL, NPl , NMl ,N2 ,NUl ,NLi , IJ , IK 
COMMON /GEOM/ X E ( 66 ) , XM { 66 ) , XC ( 66 ) , Z UB ( 5 1 ) ,ZLB I 51 ) , Z A ( 66 ) , Z B { 66 ) , 
1 ZEUI66) ,ZEL (66 ) , ZZI66 ) , ZUC31 ) , ZL (51 ) 

COMMON /PAN/ WUB(5GJ ,WL0(50) , WC3 ( 5 0 ) , WU ( 66 ) ,hL(66) ,B£TA,ALF8 
COMMON /PAN2/ CPU ( 66 ) , CPU66) , UU(66) ,UL( 66 ) , I TERN 
COMMON /SLR/ HU( 50 ) ,HL (50) ,THETAU( 50 J , THET AL (50) ,DELSU( 50 ) , 

1 DELSL(5G),SEPX,HSEP,CPIN'/(50),NJ 

DIMENSION TEMPDU(50), T£MP0L(5G) 

ALFAR= ALFA/57.29578 
CALL GCNV 
THETAU ( 1 ) = 0 ,00 5 
THETALI 1 )=0.005 
HU ( 1 ) = 1 . 5 
HL( 1)=1.5 

UIi\F=0. 135«1C85. 104 
1= 1 

5 1=1+1 

IF ( I .GT ,NA ) GO TO 2 5 
hUB = HU( I-l ) 

TH£TAtt=THETAU( I- 1 ) 

F=0.025>!=HU ( I-l )-U.022 

CF = 0 .058*( ( C.9 3- 1. SS’i'ALOGlOI HU ( I-l ) ) ) «« 1 . 705 ) * ( ( 2200 0 . «T FET AU ( I 
1 -1 ) l^^I-O.ZSS) ) 

OOOX=J. 5«CF- ( THETAU ( I-l I^i'CJUI I )-UU( I-l ) ) * ( HU ( I- 1 ) + 2 . ) ) / ( UU ( I ) 

1 '!'(XC( I )-XC( I-l) ) ) 

CHOX=( 0 .^^l-HU ( I-l )«(HU( I-l )-!.)/ THETAU ( I-l) ) - ( CF- ( F* ( HUI I-D-l. ) 

1 /HU( I-l ) ) - ( 2 .-THETAUI I-l )* (UU( I )-UU( I-l ) )=i' (HU( I-l. ) + l. ) 

2 / ( UU ( I ) -x (XC ( l)-XC( I-l) ) ) ) ) 

THETAUI I )=THETAB+DGOX«(XC( I )-XC( I-l) ) 

IF ( THETAUI I ) .LE .0.0) THET AU ( I ) = TH E T AU ( I - 1 ) / 2 . 0 
HU( I ) = HUB + DHCX=!=( XC( I )-XC ( I-l) ) 

PRINT H-VALUE 

WRITE(6,8) HU( I ) 

3 FORMAT! IX ,F10 .5 ) 

IF IHU(2).GT,I,6 ) HU(2)=1.45 
IF (HU( I ).LE.l .0) HU(I)=1.05 
IF (hU(l),GE.3.0) HU(l)=HSEP-.3 
IF ( XC( I ) .LE.15) GO TO 5 
IF ( HU( I ) .GT.HSEP) GO TO 10 
IF (I .Eg. I SEP) GO TU 15 
GO TO 3 
C 

10 IF (I.LT.ISEP) ISEP=ISEP-1 
I SEPT= I 
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CPUP=CPU( I SEP+l ) 

GO TO 20 

15 IStP=lSEP+l 

CPU( ISEP)=2.-CPU( ISEP-l)-CPU(ISEP-2) 

CPU( IS£P)= (CPU( ISEP) +CPUP)/2. 

ISEPT=I 

20 SEPX = XC( I SEP) 

25 J=2 

IF (CPL(J-l) .GT. 1.0) CPL(J-L) = 1.0 
!30 J = J + 1 

IF ( CPL ( J- I ) .GT . 1 .0 ) CPLIJ-l)=1.0 

IF { ICPLI J-l )-CPL IJ-2) ) .GE.0.0 ) GO TO 30 

NJ = J 

PRINT LOCATION OF THE AIRFOIL FRONT STAGNATION POINT 

aR I IE( 6 ,33 ) i\J 
33 FORMAT! lOX, I 3 ) 


HL ( NJ- 1 ) =1 .4 
THETAL (NJ-1 )=0 .005 
C 

OJ 35 I=NJ,NL 
THETA9 = THETAL( I-L) 
hLB = HL ( I-l ) 

F = 0.025^-HL ( I-l )-C.022 

CF = 0 .OSS’:* ( ( 0.93-1. 95 oALOGl 0 ( HL I I- I ) ) ) . 705 ) “J* ( ( 22 000. ^THETAL ( I 

1 -I ) )’i‘«(-0.268) ) 

0O0X=0 .5*0 F- (THETAL ( I- 1) ^lOLII )-UL {I-l))*(HL(I-l) + 2.))/(CL(I) 

1 *(XC(I )-XC( I-l ) ) ) 

CHOX=( 0 .5*HL ( I-l )*(HL( 1-1)-1.)/THETAL( I-l) ) * ( CF- ( F* ( HL I I-D-l. ) 

1 /HL( I-l ) )-(2 .^i'THETAL! I-l )*(UU I )-UL( I- I ) ) * ( HL ( I- 1 . H- 1 . ) 

2 /( UL(I ) *(XC! I )-XC( I-l) ) ) ) ) 

THETAL ( I ) = ThETAB+OOOX*{XC! I )-XC( I- 1 ) ) 

IF ( THETALd ).LE.0,0) TH ET AL ( I ) = TH E T AL ( I- 1 )/ 2 . 0 
FL ( I )^HL6-t-DH0X* ( XC ( I )- XC ( I- 1 ) ) 

IF ( HL ( i\J) .GT. I .6 ) HL(NJ)=1.45 
IF ( HL ( I ) .LE.l .0 ) HL(I) = 1. 05 
IF ( HL ( I ) .GE .3 .0 ) HL ( I ) =HSEP-0.3 
35 CONTINUE 
C 

THET AL I NL )=THETAL( NL-i ) + ( THET AL ( NL- 1 )- THE T AL ( ML- 2 ) ) *( XC ( NL) -XC (ML 
1 -1) )/(XC{NL-l)-XC(NL-2)) 

HL (NL) = hL( NL-1 ) MhL (NL-1 )-HL (NL- 2 ) ) * ( XC ( NL )- XC ( N L- 1 ) ) / ( XC ( N L- 1 ) 

1 -XC(NL-2)J 

C 

00 40 1=1, ISEPT 
40 OELSCT I )=HJ( I) * THETAU( I ) 

C 

ISEPT1=ISEPT+L 
00 .4 5 I = ISEPTi ,N 
45 OELSU I )=OELSU( ISEPT) 

C 

DO 5 0 I=.\'J,ML 
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c 


50 DELOLI I ) = HL { I ) *THETAL( I ) 
^SJ=^A/2 


5o 

C 

60 

C 


65 

C 


70 

C 

75 

80 

C 

85 

C 


90 

C 


C 

95 

C 


C 

100 

C 

105 

C 

lie 

c 

115 

C 

120 

C 


OJ 55 I=1,NSU 

IF (DELSU( I) .GT.CELSU( I + l) ) DELSU( I )=OEL SU (I +1 )• 
\'SL= {.■■U + IMJ ) /2 


DO 60 l=iNIJ,NSL 

IF ( DElSLI I ) .OT.OELSLI 1+ 1) ) DE L S L ( I) =0 E L SL ( H-1 ) 
IF ( I lEWN.LF.A) GO TQ 75 

IF ( ISEP.GT. ISEPG) TEMPDUI I SEP ) =DE LSU( I SEP ) 
on 6 5 1 = 1,1 SEP 

OELSOI I ) = ( OELSLK I ) +TEMPDUI 1)1/2. 

IF ( NJ .LT .NJG) TE;^PDL(iNJ)=DELSL( :NJ) 

DO 70 I=NJ,)\L 

DELSLI I ) = ( DELSK I M-TEMPOU I ) )/2. 

CO 80 1 = 1 , ISEPT 
Z£U( I )=ZU( I ) +DELSU( I ) 

00 35 I = I SEP T1 , I SEP 
ZEU( I )=ZU{ I ) +DELSU( ISEPT) 

ISEP1=ISEP+1 

CO 90 I=ISEPl,N 

ZEU( I)=ZEU{I)fOELSU(ISEPT) 

DELSLI i\J ) = 0. 005 
0ELSL(i\J+l)=0.01 
0ELSL(N.J + 2) = C.C15 


DO 9 5 l=i\J,i\L 

ZEL ( I ) =ZL I I )-DELSL I I ) 

NL 1=\L+1 

IF Iin.LE.NL) go to 105 

CO iOO I = NLl,i'i 

ZELI I )=ZEL I I )-0ELSL ( \L ) 

(M0=ISEP 
NU1 = NU-H 

DU no i=.nj,NA 

CPU! I ) =CPUI I SEP) 

CO 115 I=1,'JU 
lEMPDlJI I )=DELSU( I) 


DO 120 nf'lJ.NL 
TEMPCL ( I ) = DELSL ( I ) 

IoEPG= I SEP 
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N JG=NJ 


CALCULATE .m6w SLCPES 
CALL AFSL 

WL (9) = ( WL( 10 )+WL (d ) )/2. 

S^^.UUTH CUT LAST TEN PANEL SLCPES 

CALL LEASO 
ITERN=ITERN+l 

RETURN 

END 
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SUdRCUTIKF CCNV 

FINC THt URDIi'NiATtS AT THE CCNTRCL PCINTS 

CG.TMCN /REAO/ ATI rLE(20) ,N , N A , N T ,X y I , U M T , R , XT E , XO , X F , X C PT , WT E 1 , 

1 XA(‘31 ),Zr(51 ) ,ZC{5l),ELI66) ,ALFA,VvUENO 

COMMON /IN T/ 1 U , IL , I S , ILAM, I SEP ,NU ,NL, NPl ,NMl ,N2 ,NU1 ,NL1 , IJ , IK 
COMMCN /GEOM/ XE{66)»XMI66),XCl66),ZUB(51),ZL3(5I),ZA(66),Zb(66), 
1 ZEUI66) ,ZEL(b6) ,ZZ(66),ZU131)»ZL(51) 

COMMON /PAN/ WLliil 50) ,WLB{ 50 ,WCB( SO) ,WUI 56) ,V^L (66) ,BETA, ALFB 
common /PAN2/ CP0(66 ) , CPL(66 ) , UU( 65) ,UL( 66) , I TERN 
C 

DO 5 J=1,N 
5 ZZ(J)=ZEO(J) 

C 

MUCH=N-2 
K=1 
iO 0=1 
1=0 

15 1=1+1 

IF (I.GT.MUCH) GO TO 35 
C 

20 TERM1={ZZ( I+1)-ZZ( I ) )/(XE( I+1)-XE( I) ) 

TERM2={ZZ{I+2)-ZZ(I+l))/(XE{I+2)-XE(I+l)) 

C= ( TERM2-TEPM1 ) / (XE ( 1+2 )-XE( I ) ) 

B=TERn-C«(XE( I ) +XE( I+l) ) 

A = ZZ ( I )-d-+-XE ( I )-C*X£ ( I ) *XE( 1 ) 

C 

25 XI=XC(L) 

IF ( XI .LT .X£ ( I ) ) GO TC 30 
IF (XI .EQ.XE ( I ) ) GC TO 50 
IF (XI .GF.XE ( I + l )) GO TO 15 
GO TC 35 
C 

30 1=1-1 

IF ( I .GT.O) GO rc 20 
Z I = 0 . 0 
GO TC 40 
C 

35 ZI=A+f3«XI+C<^XiyxI 
40 IF (K.EG.l) GO TC 65 
IF (K.EG.2) GO TO 7C 
C 

45 IF (O.GE..M) GO TO 55 
L = L+1 
GO TC 25 
50 ZI = ZZ( I ) 

GO TO 40 

55 IF (K.EQ.2) GO TC 75 
K=K+ 1 
C 

DO 60 J=I,N 
60 ZZ ( J ) =ZEL( J ) 

GO TO 10 
C 

65 ZEU( L) =Z I 
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GO TC 45 
70 £EL( L) =ZI 
GO TO 45 

7 5 IF ( I go to 30 

GO TO 90 

80 DO 35 1=1, NA 
ZU( I ) = ZEU( I ) 

ZL ( I ) = ZEL( I ) 

85 CONTINUE 

90 SETU3N 
END 


TOO 
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SUbRGU r INE LEASG 

SMCUTH CLT THE LAST TEN F*Ai\EL iLGPES ON THE LOWER SURFACE 

CCM^<0,^1 / INT/ 10, IL ,1 S , ILAN , 1SEP,NU,NL, NPI,N:U tN2 ,NLl »NL 1 ,1 J , IK 
COMMUN /GEOM/ XE ( 66 ) , X,M ( 6o ) , XC ( 66 ) , ZUB ( 5 1 ) , Z LB ( 5 1 ) , Z A { 6 6 ) , Z ri ( 66 ) , 
1 Zt0166) ,ZEL166) ,2Z(6 q) ,ZU131) ,ZU5U 

COMMON /PAN/ a/UB (50) , WLB ( 50) ,WCB (50 ) ,WU (66 ) , aL ( 66) ,BETA, ALFB 
COMMON /3LR/ HO ( 50 ), HL ( 5 0 ) ,THETA J{ 50 ) , THE T AL ( 5 0 ) , DE L S J ( 5 0 ) , 

1 0ELSL(50),SEPX,HSEP,CPI.NV(50),NJ 

OIMENSIGN AA(30,30), 66(30,30), CC(30,30), BC(30), TL(30), X(30), 
1 Y(30) 

NN = NL 
NS = NL-0 
M=2 
L = l 

DO 5 I = NS,:MN 
X ( I) =XC ( I )-XC( NS ) 

5 Y{ I )=>,L( I )-WL( NS) 

DO 10 I=NS,NN 
10 CC( I ,I )=X( I ) 

DO 15 0 = 2, M 
CO 15 I=NS,NN 
15 CC( I ,J)=CC( I ,0-1 )>!‘X( 1) 

CO 20 0=1, M 
00 2C I=1,M 
AA ( I , J )=0.0 
DO 20 K=NS,NN 

20 AA( I ,0 ) =AA( I ,0 ) +CC ( K, I )#CC( K ,J ) 

DO 2 5 0 = 1 ,L 
DO 25 l=l,M 
BB( I , J )=0.0 
DO 2 5 K = NS,N,\ 

25 3B ( I , 0 ) = BB ( I ,0 ) +CC ( K , I )*Y( K ) 

SOLVE FCR SLOPE COEFFICIENTS 

CALL SIMEQ ( A A , t3 ti , HC , M , 3 0 ) 

DO 3 5 0 = NS , Nl\ 

TL( 1 )=X(0) 

SOM=au( 1 )*TL ( 1 ) 

DO 3C I=2,M 
TL ( I )= TL { I-l )'-i=X ( 0 ) 

SUMT=BC( 1 ) -TL( I ) 

30 SUM=SUM+SUMT 
35 Y( 0 ) =SLM 

CO AC 1='\S,NN 
AO WL ( I ) =V ( I ) ^-,■^L( i\S 1 

P£ Tu;^, 

END 
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SUdRCUFIinE AFSL 

CALCULATE FhE SLCPE OF THE AIRFOIL FROM ITS ORDINATES 


C 


C 


C 


C 


C 


C 


C 


C 


c 


COMMON /READ/ A T I F L E ( 2 C ) , N , N A , \T ,X M I , UMT , R , XT E , XO , X F , X C PT , WT E I , . 
I XA ( 31 ) , ZT (31 ) , ZC( 51 ) , EL (66) , aLFA , WUEND 

COMMON /INF/ I J , IL , I S , ILAM, ISEFtNU ,NL, NP 1 , NM 1 , N2 ♦ NU 1 1 NL 1 , IJ , IK 
COMMON /GEUM/ XE ( 66 ) , XM ( 66 ) , XC ( 66 ) t ZUB I 5 1 ) , ZLD I 5 1 ) , Z A { 66 ) , Z 8 ( 66 ) , 
1 ZEU(66) ,ZEL (66 ) ,ZZ(66) , ZU(31 ) , ZL I 51 ) 

COMMON /PAN/ WUB ( 50) i W L6 ( 50 ) 08 ( 5 0 ) » WU ( 66 ) , NL ( 6 6 ) , l3 ET A , ALF 3 

COMMCN /PAN2/ C PU ( 66 ) , CP L ( 66 ) , UU ( 6 6 ) , UL ( 66 ) , I T E R N 

COMMON /BLR/ HU ( 50 ) , HL ( 50 ) , THE T AU ( 50 ) , THET AL ( 50 ) , 0 EL SU { 50 ) , 

1 DELSLI50) ,SEPX,HSEP,CPINV(50) ,NJ 

DIMENSION DZDX(5G) 


ALPHR = ALFA«3. 16159/ 130 . 
□0 5 J=1,NU 

5 ZZ(J)=ZEU(J) 

« 

K=1 

MUCH=NU-2 
1 = 0 

10 1=1+1 

IF ( I.GT.MUCH) GC TO 15 


TERM1= UU I + 1)-ZZ( 1 ) ) / ( XC( I + 1)-XC( I ) ) 
FERM2=(ZZ( I+2)-ZZ(I+l)J/(XC(I+2)-XC(I+l)) 
C= ( TERM2-TERM1 ) / ( XC ( I +2 )-XC ( I ) ) 
c 3=TERM1-C« (XC( I ) +XC( I + l ) ) 

A=ZZ(I )-o«XC(I )-C^''XC(I )*XC(I) 

15 ZI=A + d«XC( I ) +C*XC( I )*XC( I) 

OZnx( I )=B + 2.«C=!‘XC( I ) 

IF (K.NE.l) GO FC 20 
IF ( 1 .CE.NU) GO TO 25 
GU FC 10 

20 IF (I .GE.NL) GC TO 40 
GO TO 10 

25 IF ( K. EG.2 ) GO TC 40 


DO 30 I=2,NU 
30 'VJ( I )=OZOX( I ) 


K = K+ 1 

GO 3 5 J = NJ » NL 
35 ZZ(J)=ZEL(J) 


MUCH=NL-2 
I=NJ-1 
GO TO 10 


40 DO 43 I=NJ,NL 
45 WL( I )=OZGX( I ) 


RE FUR\ 
EN 0 
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SUdR JU r INE P AR INB ( X , Y , I 1 , 1 2 , X 0 , X I N , YGUT » OY OUT , Y lOU T » I XXX) 

THIS SUBROUTINE INTERPCLAFES S E T WE EN PC I N F S BY PASSING A PARABOLA 
THROUGH THREE AJJACENT POINTS WITH THE DESIRED POINT IN THE SECGN 
INTERVAL UR CUINCIDING WITH THE MIDPOINT. INTEGRATION OCCURS OVE 
THE INTERVAL BETWEEN XO AND XIN, BOTH XO ANC XIN MUST LIE WITHIN 
the TABLE VALUES X(Il) AND XII2). 

DIMENSION XI 1) , Yd) 

C 

1 = I XXX 
Y0UT=0.0 
DY OUT =0.0 
YIOUT=0.0 

IF ( II .EG. 12 ) GO TO A5 
IF { I2-I 1 .EQ.l ) GO TO 50 
IF {XII U.GT.XI li+L) ) GO TO 10 
IF (XII )-XIN) 5,20,20 
5 CONTI -JOE 
GO TO 20 
C 

10 DO 15 1=11,12 

IF (X( I )-X IN ) 20 ,20, 15 
15 CONTINUE 
C 

20 I=I 

IF { XIN.NE.XI I ) ) 1 = 1-1 
IF ( I .LE . II ) I = I 1+1 
IF II.GE.I2) 1=12-1 

IF { Y( I +1 ) .EG. y ( I) .OR. Y( I ) .EQ. Y( I-l) ) GO TO 55 
IF (Y( t+1) .EG.YI I-l) ) GO TO 55 
C 

A= (X ( I +1 )-X( 1 ) ) / ( Y( I + l )-Y( I ) )-( XI I )-X( I- 1) ) / {Y( I )-Y ( I-l ) ) 

A=A/( Y( I+l )-Y{ I-l ) ) 

B=(X( I )-X( I-l) ) / IY( I )-YU-l) )-A=:«(Y( I ) + Yl I-l) ) 

C = X ( I )-0*Y ( I )-A<=Y( I ) *Y ( I ) 

IF {A. EG. 0.0) GO TC 55 
C 

D=-B/ ( 2 .G<=A) 

IF I 3 . E G . 0 . 0 ) GO T C 2 5 
E = 4. 0-A/ ( B’i'd ) 

F=1 .0-E*IC-X IN ) 

G= 1, J-E=!=l C-XC) 

IF ( G.LT.O .O.UR.F.LT.O.O) GO TO 55 
C 

Yl = D-->( 1.0-SUKT (F ) ) 

Y2=D+-( 1 .0 + S0KT(r ) ) 

IF ( Yl I ) .L T. Y( I- 1) ) GO TC 40 

IF ( Yl .GE. Y I I-l ) .AND.Y l.LE.Y ( I +1) ) GO TO 30 
IF ( Y2.GE. Yd-1 ) .AND.Y2.LE.Y( I +l ) ) GO TO 35 
GO TO 55 
C 

2 5 D=2.C/( B,0#A ) 

E=XIN-C 
F = X C-C 

IF { E.LT.O.J. OR. F.LT.Q.O) GO TO 55 
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YGUT=S:RT( E/A) 

DYUUT= 1. J/ ( 2.0=!'A’:^YUUT) 

YI Ob T = G=> { SQRT( E*E-"t )-SQRT(F’«=F=:=F) ) 

GO r G o o 

c 

30 YUUT=Y1 

DY0bT= L.0/(2.0*A*YH-3) 

YI0UT = D*(X Ii'j-X0-2.0=)‘( SQRT( F*F«F ) -SQR T I G«G=i'G ) ) / ( B.O’i'E ) ) 
GO TC 60 
C 

35 YUUT^Y2 

OYOU r= 1.0/ { 2. J«A*Y2 + B ) 

YI0UT=D« (XIi\-X0 + 2.0«{ SORT( F#F*F )-SQR I( G*G*G) ) / I 3.0^E ) ) 
GO TC 60 
C 

40 IF ( Y1 ,LE.Y( I-l J .AiMU.Yl.GE.YI I + l ) ) GO TU 30 
IF {Y2.LE. Y( I-l) .AN0.Y2.GE.Y(H-1) ) GO TO 35 
GO TC 55 
C 

45 DYObT=J.C 
YOUT=Y (ID 

YIQUT = Y( I 1 )^= (X IN-XO) 

GO TO 6C 
C 

5 0 OYCUT=( Y( I 2)-Y( I 1) ) /( X{ I2)-X( I 1) ) 

YUUT = Y ( I 1 )+DYaUT>"^ ( X IN-XI I 1 ) ) 

Y0 = Y( I 1 ) ^-OYObT’f'I XG-X( ID) 

YlOUT=0.5«( YO+YOLT)=i'( XIN-XO) 

GO TO 60 
C 

55 CALL PARI.NT ( X , Y , 1 1 , I 2 t XC r X I K , YOOT ,OYCUT , Y I ObT ) 

60 RETURN 
C 

END 
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SUBRCUT I mE P AR I N r ( X »Y , I L , 12 ,XOf XI , YCur^DYOUr , Y IGUT ) 

THIS SUbRUUriNE INTERPGLATES 8 ETaEEiN POINTS 8Y PASSING A PARABOLA 
THROUGH THREE AOJACEN'T PCI, MTS !«irH THE DES I R hO PC I NT IN THE SECON 
INTERVAL OR CJI.JCIOING WITH THE MIDPOINT. INTEGRATION uCCORS IN 
THE I.NTERVAL BET'/^EEN XU AND Xl.N. 

DIMENS ICi\ xm , YU) 

IF ( II .ED. 12 ) GO TO 55 
IF ( I2-I I .EQ .1 ) GO IQ 50 
IF {X( ID.GT.XI IU-1) ) GO TO 10 

DO 5 1=11,12 
IF (X( n-XIN ) 5, 20,25 
5 CONTINUE 
GO TO 60 
C 

10 DO 15 1=11,12 

IF (Xf I )-X IN ) 25,20, 15 
15 CONTINUE 
GO TO 60 
C 

20 I = I 

IF ( I .EG. I 2) GO TO 25 
GO TC 30 
C 

25 1=1-1 

30 IF (I-Il) 50,35,40 

35 A=(Y(I+2)-Y(I+l))/{(X(l+2)-X{I+l))«(X(If2)-X(I)))-{Y(I+l)-Y(I)) 

1 /( ( X( I-t-i )-X( I ) ) «(X I 1*2 )-X( I ) ) ) 

GO TC 45 
C 

4 0 A= { Y( I *1 )-Y{ I) ) / ( (XI I *1 )-X( I > ) - (X( I + D-X ( I-l ) ) )-( Y( I )-Y ( I-l ) ) 

1 /((Xm-X{I-l))*lX(I + l)-X(I-U)) 

45 o=(Y(I+l)-Y(I) )/(X(I+l)-X(I))-A«(X(I*l)+X(I)) 
c=Y( I )-A«x( I )*x( n-a*x{ I ) 

c 

YOUT=A*XIN«X IN+3’f=XIN*C 
DYOUT=2.C-A«XIN*B 

YIOUT=(X I,',=!=-3-X0>i‘-3 ) *A/3.0 + (XIN^*2-XOX^«^2 )’^3/2.0 + C* ( XIN-XC) 

GO rC 65 
C 

5C UYOUT=( Y( I l+D-Y n 1 ) I / I X ( 1 1 + 1 ) -X ( I 1 ) ) 
yUUr = Y ( I 1 ) +DYOUTV 1 XIN-XI I 1 J ) 

Y0=Y ( I 1 ) +DYUUT* { XO-X (ID) 

YIOLT= ( YUOT + YO ) * ( X I N- XO ) /2 . 0 
GO TO 65 
C 

55 DYGUT=0.0 
YUlJT=Y (11) 

YIOUT=Y{ ID- ( XIN-XO) 

GO TC 65 
C 

6 0 DYGLiT=lY( I 2 ) - Y ( I 2- 1) ) / { X ( I 2 ) -X ( I 2- 1) ) 

YuJT=Y{ I2D-DYOUT-(XI;V-X{ 12) ) 

VO=Y ( I 2 )+UYCUr=(= ( Xj-X (12)) 

Y I GO T = ( YC'AJ T + YO ) X I N- X C ) / 2 . 0 
C 

65 RETURN 105 

END 



non 


SUB^CUTI.NE SIMEQ ( A , Y , X , , NM AX ) 

THIS IS THE STANCARC I EM ROUTINE TO SOLVE SIMULTANEOUS 
DIMENSION A(NMAX,1), Y(l), XII) 

M = N-l 

DU L5 I=1,M 
AI I=A( I , II 
L=H-l 

DO 15 J=L,N 
AJ I = A( J , I ) 

IF (AJI ) 5,15,5 
5 DO 10 K=L,N 

10 A( J,K)=A( J,X)-A( I,K)*AJI/AI I 
Y( J)=Y( JI-Y( n«AJI/AI I 
15 CONTINUE 

X(N)=YIN)/A(N,N) 

DO 25 r=l,M 
K = N- I 
L = K+ 1 

DO 20 J=L,N 

20 Y( K) =Y ( K )-X ( J) ^=A(K,J) 

25 X(KI=Y(K)/A(K,X) 

RETURN 

END 


QUA nONS 


; 


f 


106 



